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CHAPTER I 


INTRODUCTION 

Structural failure is rarely a "sudden death" type of 
event, such sudden failures may occur only under abnormal 
loadings like bomb or gas explosions and very strong 
earthquakes. In most cases, structures fail due to damage 
accumulated under normal loadings such as wind loads, dead and 
live loads. The consequence of cumulative damage will affect 
the reliability of surviving components and finally causes 
collapse of the system. The cumulative damage effects on 
system reliability under time-invariant loadings are of 
practical interest in structural design and therefore will be 
investigated in this study. 

The scope of this study is, however, restricted to the 
consideration of damage accumulation as the increase in the 
number of failed components due to the violation of their 
strength limits. Progressive failure processes such as 
corrosion, fatigue and crack growth are not investigated in 
this study. 


1 
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1.1 Background and Significance of the Study 

Structural designs have been traditionally based on 
deterministic design methodology. The deterministic method 
considers all design parameters to be known with certainty. 
This methodology is, therefore, inadequate to design complex 
structures subjected to a variety of complex, severe loading 
conditions. These complex conditions introduce uncertainties 
and so the actual factor of safety remains unknown. In the 
deterministic methodology the contingency of failure is 
discounted, so there is a use of a high factor of safety. 

Probabilistic design method is concerned with the 
probability of non-failure performance of structures or 
machine elements. Probabilistic methodology is a convenient 
tool to describe, or model, physical phenomena too complex to 
treat with the present level of scientific knowledge. It is 
much more useful in situations where the design is 
characterized by complex geometry, possibility of catastrophic 
failure or sensitive loads and material properties. 

1.1.1 Comparison between Deterministic and Probabilistic 

Design Methodology 

The probabilistic design methodology produces 
designs that are robust and allows the quantification of the 
level of reliability in the design, as opposed to 
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deterministic designs. Hence, it is begining to attract more 
attention than the traditional deterministic design. 

Probabilistic design procedures promise to improve the 
cfuslity of engineered systems for the following reasons: 

1 . Probabilistic design incorporates given statistical 
data explicitly into design algorithms. Conventional design 
discards such data. 

2. It is more meaningful to say, "This system has a 
probability of 10~ 4 of failing after 1000 hours of operation, " 
than to say, "This system has a factor of safety of 2.3." 

3. Rational comparisons can be made between two or more 
competing designs for a proposed system. Without other 
considerations, the engineer chooses the design having the 
lowest probability of failure, or basis for developing 
economic strategies. 

4. An "optimal" design of a system results when each 
component chosen so that its probability of failure is the 
same . 

5. By treating each nonstatistical uncertainty as a 
random variable, its effect on the final design can be 
quantified. 

6. Probabilistic-based information on mechanical and 
structural performance can be used to develop rational 
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policies toward pricing, warranties, etc. 

1.1.2 Structural Reliability under Time-invariant Loads 

This study primarily focusses on the effects of time- 
invariant loads on the structure. The effects of time- 
invariant loads on element and system reliability are 
discussed below. 

1.1. 2.1 Element Reliability 

The study of element reliability under cumulative damage 
is to include the system effects into element reliability. In 
the current codes such as CEB[1], LRFD [ 2 ] and AASHT0[3] 
specifications, the design of a structural system goes through 
the design of components and connections individually. The 
target element reliability and safety are achieved by making 
them satisfy the limit state functions of local strength with 
a high degree of probability. What is the reliability of the 
individual component once it is in the actual configuration? 
How do the system effects influence the element reliability 
and which components are more vulnerable than the others? What 
impact do these questions have on a current reliability-based 
design code, like the AISC load and resistance factor design 
(LRFD) code? These are some preliminary questions sought to be 
answered in this study. 

Mahadevan and Haidar [4] used the Stochastic Finite 
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Element Method (SFEM) to investigate the magnitude of system 
effects on component reliability in framed structures designed 
by the LRFD code. However, their analysis was based on linear 
elastic behavior. The effect of geometric non-linearity was 
included in SFEM-based reliability analysis by Liu and Der 
Liureghian [ 5] and Haidar and Zhou[6], The effect of material 
nonlinearity has been considered by many researchers to 
estimate the overall system reliability, but the focus of this 
study is the reliability of individual elements affected by 
the formation of plastic hinges elsewhere in the structure. 
Therefore, a rational procedure has been developed in this 
study to account for the effect of structural system damage. 

Although system reliability research has been active for 
the past twenty years, it still has not been applied in 
practical design. The inclusion of system effects on element 
reliability may offer a solution to this problem so that the 
element-based design can account for system effects. 

1.1. 2. 2 System Reliability 

The collapse of a system is the culmination of cumulative 
damage of components. This idea resulted in the development of 
several failure path identification techniques including 
branch and bound method, (5-unzipping method, etc. However, 
these techniques are difficult to implement in case of large 
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structures, because they are time consuming. This is one of 
the important reasons for the slow application of system 
reliability in modern design[7]. This study and analysis are 
used to examine the performance of the LRFD approach in the 
design of realistic structures, resulting in several important 
observations . 

In this study, the loadings are idealized as time- 
invariant. In other words, the reliability so obtained 
corresponds to that under one load application though it may 
represent some extreme value of the load over a given period. 
However, the reliability of an element or a structure varies 
over its lifetime, due to repetitive load applications causing 
accumulated damage, degradation of material resistance over 
time, corrosion, wear etc. 

1.2 Research Objectives and Organization of the Report 

The above discussion leads to the following objectives: 
1. Discussion of the probabilistic design methodology 
in depth and an overview of the software, NESSUS (Numerical 
Evaluation of Stochastic Structures Under Stress) used 
primarily in this project. This is described in detail in 


Chapter II. 

2. Discussion of the LRFD source codes and their 
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application to this project. This is described in Chapter III. 

3. Development of a failure path-based procedure to 
estimate system reliability of assumed steel structures, along 
with the development of a computational procedure to estimate 
the element reliability under cumulative damage. This is 
described in Chapter IV. 

4. The results of the system reliability analysis, their 
interpretation and explanation, are described in Chapter V. 

5. The summary and conclusions of the present study and 
suggestions for further research are presented in Chapter VI . 

6. The appendix A lists the computer program formulated 
for lognormal distribution. Appendix B gives the detailed 
loading calculations done for the structures in accordance 
with the Uniform Building Codes. Moment analysis of the 
structures, which is done by finite element software (STAAD- 
III) is given in appendix C. The algorithm and flowchart to 

operate NESSUS for probabilistic design is listed in Appendix 

D. 



CHAPTER II 


PROBABILISTIC DESIGN METHODOLOGY 

2.1 Role of Probability in Engineering 

Quantitative methods of modeling, analysis, and 
evaluation are the tools of modern engineering. Some of these 
methods have become quite elaborate and include sophisticated 
mathematical modeling and analysis, computer simulation, and 
optimization techniques. However, irrespective of the 
sophistication in the models, including experimental 
laboratory models, they are predicated on idealized 
assumptions or conditions; therefore, information derived from 
these quantitative models may or may not reflect reality 
closely. 

In engineering designs, decisions are often required 
irrespective of the state of completeness and quality of 
information, and thus are made under conditions of 
uncertainty. In other words the consequence of a given 
decision cannot be determined with complete confidence. 
Besides the fact that the information must often be inferred 
from similar circumstances or derived through modeling. Many 
problems in engineering involve natural processes and 
phenomena that are inherently random; the states of such 
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phenomena are naturally indeterminate and thus cannot be 
described with definiteness. For these reasons, decisions 
required while engineering planning and design invariably must 
be made, and are made, under conditions of uncertainty. 

The effects of such uncertainty in design and planning 
are important. To be sure, the quantification of such 
uncertainty and evaluation of its effects on the performance 
and design of an engineering system, should include concepts 
and methods of probability. Further, under conditions of 
uncertainty, the design and planning of engineering systems 
involve risks, and the formulation of related decisions 
requires them to be risk free. The problems of uncertainty in 
the design can be overcome by applying the methods of 
probability. Thus, the role of probability is quite pervasive 
in engineering. It ranges from the description of information 
to the development of bases for design and decision making [8]. 
Many phenomena or processes of concern to engineers contain 
randomness, that is, the actual outcomes are sometimes 
unpredictable. Such phenomena are characterized by 
experimental observations that are different from one 
experiment to another, even if performed under identical 
conditions. In other words, there is usually a range of 
measured or observed values and within this range certain 
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values may occur more frequently than others. Clearly, if 
recorded data is of a variable exhibit scatter or dispersion, 
the value of the variable cannot be predicted with certainty. 
Such a variable is known as random variable and its value or 
range of values can be predicted only with an associated 
probability. When two or more random variables are involved, 
the characteristics of one variable may depend on the other. 

Since there is a range of possible values of random 
variable, we would be interested in some central value, such 
as the average. In particular, because the different values of 
the random variable are associated with different 
probabilities, the weighted average is taken into 
consideration. This weighted average is known as sample mean 
value of the random variable. Therefore, if X is a discrete 
random variable, then the mean value, p x is obtained as 
follows 


where. 



p x is the mean 

X is the random variable. 

n is the number of observations. 


(2-1) 
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Besides the sample mean, the next most important quantity 
of a random variable is its measure of dispersion or 
variability, that is, the quantity that gives a measure of how 
widely the values of the variate are spread around its mean 
value. This deviation can be above or below its central value. 
If the deviations are taken with respect to its mean value, 
then a suitable average measure of dispersion is called the 


Var(X) 


2(* - H x ) 2 
n - 1 


( 2 - 2 ) 


variance and is computed using the following relation: 
where, 

Var(X) is the variance of the random variable X. 
Dimensionally, a more convenient measure of dispersion is 
the square root of the variance, or the standard deviation, 

o x = JVaf(Xj (2-3) 


where, 

o x is the standard deviation of the random variable X. 

Saying whether the dispersion is large or small is 
difficult, from the variance or standard deviation. For this 
purpose, the measure of dispersion about the central value is 
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more useful. In other words, the dispersion is large or small 
is meaningful only about the central value. Therefore, 
coefficient of variation (COV) is often preferred, which is a 
convenient non-dimensional measure of dispersion or 
variability. The coefficient of variation is related to the 
mean and standard deviation as follows, 


where. 


o 

COV = -i 

Vx 


(2-4) 


o x = Standard deviation of the variable X. 
U x = Mean value of the variable X. 


The application of probability is not limited to the 
description of experimental data, or the evaluation of the 
statistics such as the mean and standard deviation. In fact, 
the more significant role of probability concepts is in the 
use of this information in the formulation of proper bases for 
the design. 


2.2 Uncertainty Associated with Design 

Engineering uncertainty is not limited to the variability 
observed in the basic variables. First, the estimated values 
of a given variable (such as the mean) based on observational 
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data will not be error free. Second, the mathematical or 
simulation models (for example, formulas, equations, 
algorithms and laboratory models, that are often used in 
engineering analysis and designs are idealized representations 
of reality) . Consequently, predictions and calculations made 
from these models may be inaccurate (to some unknown degree) 
and thus also contain uncertainty. Human error can result from 
errors made by engineers and technicians during the design or 
operations phases. It can be reduced by improving the quality 
of a control program, but it cannot be avoided entirely. 
Usually, human error is very difficult to define. In study, 
human error will be treated as modeling error. In some cases, 
the uncertainties associated with such prediction or model 
errors may be much more significant than those associated with 
the inherent variabilities. 

All uncertainties, whether they are associated with 
inherent variability or with prediction error, may be assessed 
in statistical terms and the evaluation of their significance 
on the design can be accomplished by the concepts and the 
methods of probability. 
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2.3 Designing under Uncertainty 

If there are uncertainties in the design, the next step 
is to ask, how should designs be formulated or decisions 
affecting a design resolved? Presumably we may assume the 
worst conditions and develop conservative design on this 
basis. From the system performance and safety point of view, 
this approach may be suitable. However, the resulting design 
would be too costly because of over conservatism. On the other 
hand an inexpensive design may not ensure the desired level of 
performance and safety. Therefore the decisions should be made 
considering cost and safety of the design. The most desirable 
solution is one that is optimal, in the sense of minimum cost 
and maximum benefits. If the available information and the 
models to be evaluated contain uncertainties, the analysis 
should include the effects of such uncertainties. 

Probabilistic design is concerned with the probability of 
failure, or preferably, reliability. This methodology is most 
useful when uncertainties in material properties and loading 
conditions are considered. To apply probabilistic design 
methodologies (PDM) , all uncertainties are modeled as random 
variables, with selected distribution types, means and 
standard deviations. The primitive (random) variables that 
affect the structural behavior have to be identified. Every 
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design project demands some sequential stages of reflection 
before one can arrive at the final design goal. This is also 

the case with PDM. The various design stages of PDM are as 
follows . 

1. Problem Definition. 

2. Generating design parameters. 

3. Relating the defined problem to the design parameters. 

4. Data assembling and application of probability 

concepts . 

5. Probabilistic Analysis. 

6. Interpreting results. 

The design stages of PDM are shown in Figure 2-1. 



Figure 2-1 : Design stages in PDM [9] 
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1. Problem definition 

The first step a designer takes in solving a design 
problem is to find out the main objective of the design. After 
finding out the objective, the next step is to define in a 
precise manner the functional requirements, of the system or 
component to be designed. These functional requirements should 
be able to completely characterize the design objective by 
defining it in terms of specific needs. With a clear 
understanding of what one is searching for, the designer then 
goes to the next stage. 

2 . Generating design parameters 

In order to solve the defined problem, acceptable design 
parameters must be generated that will meet the defined 
functional requirements . To generate the design parameters one 
uses an appropriate design model. The various parameters like 
loads, material properties, geometry, crack size etc. are 
taken into consideration. The design parameters to be selected 
depend on the objective of the design [9]. 

3 . Relating the defined problem to the design parameters 

After defining the design parameters the designer then 
relates the functional requirements in the functional domain 
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to the design parameters in the physical domain, to be sure 
that the objective is satisfied. If the relation is 
satisfactory, the designer goes to the next stage, if not the 
relation is redefined, so that the objective is satisfied. 

4 . Data assembling and application of probability concepts 

This stage requires assembling the essential data that 
are available on the problem with regard to the design 
parameters. If some data are unavailable then it becomes 
necessary to perform a computational simulation analysis to 
generate the missing details. Once the data has been 
assembled, the next stage is to analyze the assembled data. 
NESSUS is the computer tool used to perform the analysis. 
NESSUS has three modules known as NESSUS/PRE, NESSUS/FEM and 
NESSUS/FPI . 

NESSUS/PRE is a preprocessor, which prepares the 
statistical data needed for the probabilistic design analysis. 
It allows the user to describe the uncertainties in the 
structural design parameters. The uncertainties in these 
parameters are specified by defining the mean value, standard 
deviation and the distribution type, together with an 
appropriate form of correlation. Correlated random variables 
are then decomposed into a set of uncorrelated vectors by a 
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model analysis. 

NESSUS/FEM is a general purpose finite element code, 
which is used to perform structural analysis and evaluation 

er- 

of sensitivity due to variation in different uncorrelated 
random variables. The response surface, defined in terms of 
random variables required for probabilistic analysis in 
NESSUS/FPI, is obtained from NESSUS/PRE. NESSUS/FEM 
incorporates an efficient perturbation algorithm to compute 


the sensitivity of random variables [10]. 



Figure 2-2: Modules of NESSUS 


NESSUS/FPI is an advanced reliability module, which 
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extracts the database generated by NESSUS/FEM to develop a 
response model in terms of random variables [11 ] . In this 
module the probabilistic structural response is calculated 
from the performance model. The probability of exceeding a 
given response value is estimated by a reliability method. 
Inside the NESSUS/FPI module is a sensitivity analysis 
P^ocf^sm, which determines the most critical design parameters 
in the design. The input data for NESSUS/PRE requires 
fundamental knowledge of statistics or probability theorems. 
The expected details will include determining the mean, 
standard deviation, median, coefficient of variation, 
variances etc., associated with each random variable. The 
designer also determines the probability distribution function 
that best describes each random variable. The different 
modules of NESSUS are shown in Figure 2-2. 

5. Probabilistic Analysis 

It is at this stage of the design that the designer 
defines a limit state function. The limit state function is a 
function that defines the boundary between the safe and 
failure region. In the limit state function approach for 
structural reliability analysis, a limit state function g (X) 
is first defined. The g-function, is a function of a vector of 
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basic random variables, X-(X„ X 2; X 3 jj, ) wit h g (X) . o 

being the limit state surface that separates 


— 



Figure 2-3: Illustration of Most Probable Point 

the design space into two regions, namely, the failure 
g(<0) and the safe g(>0) regions. Geometrically, the limit 
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state equation, g(X)=0, is a n-dimensional surface that may be 
called the "failure surface." One side of the failure surface 
is the safe state, g(X)>0, whereas the other side of the 
failure surface is the failure state, g(X)<0. 

The probability of failure in the failure domain O is 
given by; 

Pf = Jo Jf, (X) dx (2-6) 

where f x (X) is the joint probability density function of X and 
Q is the failure region. The solution of this multiple 
integral is, in general, extremely complicated. Alternatively, 
a Monte Carlo solution provides a convenient but usually time 
consuming approximation. The limit state function method uses 
the Most Probable Point (MPP) search approach shown in Figure 
2-3. The Most Probable Point is the key approximation point 
for the FPI analysis, therefore, the identification of MPP is 
an important task. In general, the identification of the MPP 
can be formulated as a standard optimization problem and 
solved by proper optimization methods. 

From the Figure 2-3, as the limit state surface g(X)=0, 
moves closer to the origin, the safe region, g(X)>0, decreases 
accordingly. Therefore, the position of the failure surface 
relative to the origin of the reduced variates should 
determine the safety or reliability of the system. The 
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position of the failure surface may be represented by the 
minimum distance from the surface g(X)=0 to the origin. The 
point on the surface with minimum distance to the origin is 
the Most Probable Point (MPP) . This is usually determined by 
fitting a local tangent to g(X) and moving this tangent until 
MPP is estimated. 

In the NESSUS code MPP is defined in a transformed space 
called u-space where the u's are independent to simplify the 
probability computations. By transforming g(x) to g(u), the 
most probable point, u*, on the limit state, g(X)=0, is the 
point that defines the minimum distance from the origin to the 
limit state surface. This point is most probable (in the u- 
space) because it has maximum joint probability density on the 
limit state surface. The required minimum distance is 
determined as follows. The distance from a point u*=(u : ‘, u 2 *, 
..., u n ) on the failure surface g(u)=0 to the origin is, 

D = V W 1 + U 7 + + U n (2-7) 


where, D is the minimum distance from the point on the limit 
state surface to the origin. 

The FPI code assumes only one MPP. In general, however. 
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the possibility exists that there may exist multiple local and 
global Most Probable Points. A two MPP problem can occur for 
example, if the g-function is quadratic and the search 
algorithm may result in an oscillating (non-convergent ) 
search . 

Several approaches are available to search for the MPP. 
The search procedure depends on the forms and the number of 
the g-function ( s ) . One efficient method in use is the Advanced 
Mean Value method. This method blends the conventional mean 
value method with the advanced structural reliability analysis 
method. This method provides efficient cumulative density 
function analysis and the reliability analysis. The step wise 
AMV method can be summarized as follows [12] : 

1. Obtain the g(X) function based on perturbations about 
the mean values. 

2. Compute the cumulative density function of the 
performance function at selected points using the fast 
probability integration method. 

3. Select a number of cumulative density function values 
that cover a sufficiently wide probability range. 

4. For each cumulative density function value, identify 
the most probable point. 

Another approach considered efficient as well is the 
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Adaptive Importance Sampling Method. This method focusses on 
reducing the sampling domain in the search space after the MPP 
is identified. The Adaptive Importance Sampling method is 
generally used for system reliability analysis. 

The analytical process involved in the limit state 
approach can be illustrated by a basic structural reliability 
problem. In the problem only one load effect S, restricted by 
one resistance R, is considered. 

If one considers a case when R and S are independent, 
the limit state equation can be expressed as, 

g = R - S (2-8) 

and the probability of failure can be expressed as, 

Pf = P(R-S^O) = J_j Jf R (r) f S (s) dr ds (2-9) 

For any random variable the cumulative density 
function F(x), is given by 

F x (x) = P(X £ x) = Jf x (y)dy (2-10) 

provided that x £ y 

Therefore P f is expressed as 

Pf - P(R-SsO) = Jf r (x) f s (x)dx (2-11) 

Assuming a special case of normal random variables, for 
some distributions of R and S, it is possible to integrate the 
equation (2-11) analytically and find out the probability of 
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failure. If S and R have mean p R and p s and variance's o R 2 and 
o s - respectively, the g-function has a mean p g and variance o g 2 , 
given by 

M g = Mr -M s (2-12) 

a 8 2= °R 2 + ° s 2 (2-13) 


Therefore the probability of failure is given as. 


P f = P{R-S<. 0) = P(g< o) = <t>[— 


(2-14) 


«>[ - ilpl] = 4> (-P) (2-15) 


Which reduces to. 


where (3 is defined as the safety index. 


(2-16) 


Thus the probability of failure is given as 
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Pf = 4>(-P) 


(2-17) 


which can be written as, 


Pf = 1 - 4KP) (2-18) 

The reliability of the system is given by 

P r = 1 - P f (2-19) 

where P r is the reliability of the system. 

6. Interpretation of Results 

This is the last stage in the methodology. When the 
designer approaches this stage, one interprets the results 
obtained about the initial objective. If the results do not 
satisfy the functional requirements in the stage 1, the 
designer may adjust order to achieve the set objective. 

2.4 Probability Sensitivity Factors 

In Engineering performance analysis many sensitivity 
measures can be defined. Knowing the effect of each random 
variable in the analysis is important for the designer. The 
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sensitivity information is quantified by sensitivity factors. 
Sensitivity factors suggest which random variables are crucial 
and require special attention. 

The commonly used sensitivity in deterministic analysis 
is the performance sensitivity, dZ/dX if which measures the 
change in the performance due to the change in a design 
parameter. This concept can be extended to the probabilistic 
analysis in which a more direct sensitivity measure is the 
reliability sensitivity that measures the change in the 
probability/reliability relative to the distribution 
P^^^m^ters such as the mean and the standard deviation. 
Although not automated in the code, this analysis can be 
performed by varying the parameters. 

Another, perhaps more important, kind of probability or 
reliability sensitivity analysis is the determination of the 
relative importance of the random variables. This analysis can 
be done, for example, by repeated probabilistic analysis in 
which one random variable at a time is treated as a 
deterministic variable. The results of the analyses, for 
example, are a number of cumulative density function curves or 
reliabilities. Based on the results, the relative importance 
of the random variables can be analyzed. The standard FPI 
output includes a first order sensitivity factor that provides 
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approximate relative importance of the random variables. 



CHAPTER III 


LOAD AND RESISTANCE FACTOR DESIGN FOR STEEL 

Inherent uncertainties in structural design parameters, 
such as loads, geometry, material and sectional properties, 
and boundary conditions, are well established. However, in 
traditional design procedures, these parameters are considered 
deterministic; the uncertainty is accounted for by the use of 
safety factors. Thus, in allowable-stress design method, the 
ultimate stresses are divided by safety factors to determine 
the allowable stresses. A successful design ensures that the 
stresses caused by the nominal values of the loads do not 
exceed the allowable stresses. In the ultimate-strength, or 
plastic, design, the loads are multiplied by the load factors 
to determine the ultimate loads and the fully stressed members 
are required to resist various design combinations of these 
ultimate loads [ 13 ]. 

A more rational approach to consideration of 
stochasticity in structural parameters has resulted in the 
development of the LRFD approach during the past decade. 

3.1 General Discussion on LRFD codes 

The load and resistance factor design criterion is 
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expressed by the following general formula: 

SRn^lYiQi (3-1) 

The left side of the formula relates to the resistance 
(capacity) of the structure while the right side characterizes 
the loading acting on it. 

The resistance side of the design criterion consists of 
the product 4>R n , in which R n is the "nominal resistance," and 
(f 1 is the "resistance factor." The nominal resistance is the 
resistance computed according to a formula in a structural 
code and it is based on the nominal material and cross- 
sectional properties. The resistance factor <t>, which is always 
less than unity, together with R n reflects the uncertainties 
associated with R. The factor <t> is dimension less and R n is a 
generalized force: bending moment, axial force, or shear force 
associated with a limit state of strength and serviceability. 
Interaction equations, e.g., between axial force and bending, 
may also be used to define R n for appropriate limit states. 
The loading side of the design criterion is the sum of 
products, Yi Qi/ in which Q t is the "mean load effect," and Yi 
is the corresponding "load factor." Here Yi is dimensionless 
and Qi is a generalized force (i.e., bending moment, axial 
force or shear force) computed for the mean loads for which 
the structure is to be designed. The y-factors reflect 
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potential overloads and uncertainties inherent in the 
calculation of the load effects. The summation sign in the 
equation denotes the combination of load effects from 
different load sources [13]. 

The LRFD codes were developed, based on first order 
probabilistic design methods. In LRFD, the nominal resistance 
always relates to a specific "limit state." Two classes of 
limit states are pertinent to structural design: the "maximum 
strength" (or "ultimate") limit state, and the "serviceability" 
limit state. Violation of a strength limit state implies 
"failure" in the sense that a clearly defined limit of 
structural usefulness has been exceeded, but this does not 
necessarily involve actual collapse. In case of structural 
system with "compact" beams this means that a plastic 
mechanism has formed. Serviceability limit states include 
excessive deflection, excessive vibration, and premature 
yielding or slip. 

A first order probabilistic design procedure was used to 
determine the values of 4>, R n , y and ; Q , during the 
development of the code. This is simplified method that uses 
only statistical parameters, i.e., means values and 
coefficients of variation of relevant variables and a 
relationship 3 between them, called the "safety index." 
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Probability-based LRFD criteria have been adopted in 
Canada for hot-rolled and cold-formed steel structures, the 
basic guidelines for European national codes have been 
formulated, and research on development of similar procedures 
is underway for reinforced concrete and wood structures. 
Experience gained from one effort is transmitted to newer 
projects, and the concepts of the applications of probability, 
statistics, optimization, and decision theories have become 
increasingly more sophisticated [13] . Thus the field of design 
methodology research is very active and changes occur rapidly. 

3.2 Selection of Model 

The probabilistic design format used to develop the LRFD 
criteria for steel structures is due to Cornell [14]. This 
format was selected because of its simplicity and its ability 
to treat all uncertainties in a design problem in a consistent 
manner. The format is explained briefly in the following. 

Structural safety is a function of resistance, R, of the 
structure and of the load effect, Q, acting on it; R and Q are 
random variables. An example of the definition of safety is 
given in the Figure 3-1, where the frequency distribution of 
the random variable of R-Q, called the safety margin, is shown 
and survival is indicated by R-Q, called the safety margin, is 
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shown and survival is indicated by R-Q > 0. The probability of 
failure p F of a structural element according to the 
representation of Figure 3-1 is equal to 

cr p F =P[(R-Q)<0] (3-2) 

An equivalent representation of structural safety is 
shown in figure where the probability of failure is 
p F =P[ ln(R/Q) < 0] (3-3) 

The format according to the Figure 3-1 was adopted for 
developing the LRFD criteria. 



Safety margin, R—Q Value of In (R/ Q) 

(a) PROBABILISTIC MODEL (6J DEFINITION OF SAFETY INOEX 


Figure 3-1 : Definitions of Structural Safety [12] 


If the "standardized variate” U is introduced/ in which 
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ln(R/Q)-[ln(R/Q)] m 

U = (3-4) 

°ln(R/Q) 

where, [In (R/Q) ] m and o ln (R/Q) are the mean and standard deviation 
of the natural logarithm of the ratio (R/Q) , then from 
equation 3-3, the probability of failure can be written as 
given below 

p F = P{U< -[ln(R/Q )]Ja^ Q) } 

= Fy { -[ln(R/Q)j (3-5) 

Here is the cumulative distribution function of the 
standardized variate U. The quantity [ln(R/Q) ] m /a ln(R/Q) defines 
the reliability of the element, thus it is called "safety 
index," 3- If the probability distribution of (R/Q) were 
known, 3 would directly indicate a value of the probability of 
failure. In practice, the probability distribution of R/Q is 
unknown and only the first two statistical moments of R and Q 
are estimated. In the first-order probabilistic design method 
used here, 3 is only a relative measure of reliability; a 
constant value of 3 effectively fixes the reliability as a 
constant for all similar structural elements. 


The expression for the safety index 3, 
P = [\n(R/Q)] Ja^ Q) 


(3-6) 
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can be simplified by using first-order probability theory into 


p = (3-7) 

in which and <X are the mean values of the resistance and 
the load effect, and V r and Q V are the corresponding 
coefficients of variation [ 12 ] . 

3.3 Load Combinations 

Most load effects are random functions of time. The 
following are some important load combinations to be studied: 

1. Dead load + lifetime maximum live load 

2. Dead load + sustained live load +lifetime max 

wind load 

3. Dead load + lifetime max live load + daily max 

wind load 

4. Lifetime max wind load - dead load 

5. Dead load + lifetime max snow load[13]. 

An examination of these loads follows. 

3.3.1 Live Load - Statistical information on live loads is 
usually obtained from load surveys that give the live loads in 
the particular buildings surveyed at the times the surveys 
were made. From the load combination enumerated earlier, it is 
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seen that the distribution of the lifetime maximum live load 
is also needed. Pier and Cornell [15] have modeled the live 
load as consisting of the superposition of two parts: The 
sustained live load, which remains on the floor for a 
relatively long time until an occupancy change occurs, and the 
transient live load, which occurs infrequently but with a 
relatively high intensity and short duration. The sustained 
load includes furniture and normal working personnel. The 
transient live load may be caused by people in a room. Peir 
has proposed models to derive the statistics of the lifetime 
maximum sustained load and of the transient live load. Using 
the live load models of Pier and live load survey data of 
Mitchell and Woodgate [16] , McGuire and Cornell[17] have 
derived the statistics of lifetime maximum live load. 

3.3.2 Wind Load. - There are three random variables of 
interest in case of wind loads: The daily maximum, the annual 
maximum, and the lifetime maximum wind load. Meteorological 
data are available to derive the distributions of the daily 
maximum and of the annual maximum wind speeds throughout the 
United States. The lifetime maximum wind speed is 
approximately derived as the maximum of n-identically 
distributed and statistically independent random variables 
representing the annual maximum values, where n is the 
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lifetime of the structure in years. The mean and the 
coefficient of variation of the wind load (daily maximum, 
annual maximum, or lifetime maximum) are obtained taking into 
account the uncertainties in the dynamic characteristics of 
wind and the structure. 

An analysis of 13 locations in the continental United 
States is given in the table 3-1 for a 1-yr period, which 
lists the location, the mean fastest mile daily wind speed, in 
miles per hour (V 30Qr ) , the corresponding 50-yr ANSI wind speed 
for the same location (V^) , the factor (V^/V^ ) 2 by which 
ANSI 50-yr wind pressure is multiplied to obtain the mean load 
-intensity, and the coefficient of variation of the daily wind 
speed, V TO [ 15] . 
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TABLE 3-1 Maximum Daily Wind Statistics 


Location 

V 30dm , miles 
per hour 

Vansj, miles 
per hour 

( ^ 30 Dm/^ANSI ) 2 

Vvd 

Boston 

21 

90 

0.05 

0.32 

Denver 

19 

80 

0.06 

0.38 

Minneapolis 

18 

75 

0.06 

0.33 

Chicago 

18 

80 

0.05 

0.30 

1 St. Louis 

18 

70 

0.07 

0.37 

Kansas City 

18 

70 

0.06 

0.39 

Salt Lake C 

18 

80 

0.05 

0.39 

Washington 

17 

75 

0.05 

0.36 

Dallas 

17 

70 

0.06 

0.35 

Atlanta 

17 

80 

0.04 

0.38 

Pittsburgh 

16 

70 

0.05 

0.33 

Seattle 

16 

80 

0.04 

0.37 

New York C. 

14 

80 

0.03 

0.32 


3.3.3 Load Factors- The purpose of load factors is to 
increase the loads to account for the uncertainties involved 
in estimating the magnitudes of dead or live loads. The usual 
load combinations to be considered are given below [19] . 
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1. U = 1.4D 

2. U = 1.2D + 1 . 6L + 0.5(L r or S or R) 

3. U = 1.2D + 1 . 6 (L r or S or R) + (0.5L or 0.8W) 

4. U = 1.2D + 1 . 3W + 0.5L + 0.5(L r or S or R) 

5. U = 1.2D + 1.5E + (0.5L or 0.2S) 

6. U = 0.9D - (1.3W or 1.5E) 
where, U = ultimate loads 

D = Dead loads 

L = Live loads 

W = Wind loads 

L r = roof live loads 

S = snow loads 

R = rainwater or ice load 

E = Earthquake load 

3.4 Bending Resistance of Steel Beams 

The plastic range represents the optimum capacity of the 
beam. Beams in this region are often described as compact 
beams [20]. In this range the plastic moment Mp = F y Z can be 
reached or exceeded, and this moment-level can be maintained 
for a large enough rotation so that inelastic force 
redistribution can take place and finally a mechanism can 
form. While in the elastic range of lateral-torsional buckling 



40 


the situation is clear, i.e., the member buckled or it is 
stable, the factors affecting the behavior in the plastic 
range are complex and intricately interrelated. Local flange, 
local web, and lateral-torsional distortions interact and they 
tend to build up gradually rather than form suddenly. Strain 
hardening on the one side and instability on the other side 
work against each other and they tend to balance out to give 
M=Mp at the critical length L p [21] . 

While much is known experimentally in the plastic range 
about the relationship of unbraced length and flange and web 
width-thickness ratios to rotation capacity, no generally 
satisfactory analysis that recognizes the complex 
interrelationships has yet been presented. Indeed, even if 
such a relationship did exist, its usefulness in design office 
situations would be questionable. Requiring designers to 
determine the required amount of rotation capacity to permit 
a desired level of moment redistribution would not be 
practical. The process is difficult, time consuming, and 
unreliable. Strain hardening significantly reduces the 
required rotation capacities based on ideal hinge behavior, 
i.e., = M p [ 1 9 ] . 

Studies have been made on rotation capacity requirements 
of some general structures. These studies show that for 
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practical structures, the required rotation capacity is 
small (less than two). These are usually in extreme 
structures (single-story frames with very steep gables), or in 
zones of high moment gradient, where the ideal assumptions are 
invalid. In addition, these cases usually show that at a load 
just a few percent below the maximum, the requirements are 
greatly diminished. Current rules in plastic design are not 
based on any consistent rotation capacity requirements. The 
table below shows the statistical derivations of several tests 
on beams in plastic range. 


TABLE 3-2 Statistical Data on Beam Tests in Plastic Range 


Type of member 

Number of 
tests 

(Test /predict ion) 

V p 

Statically 
determinate 
beams under 
uniform moment 

33 

1.02 

0.06 

Statically 
determinate 
beams under 
moment 
gradient 

43 

1.24 

0.10 

Statically 
indeterminate 
beams and 
simple frames 

41 

1.06 

0.07 
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3.5 Properties of Steel 

The importance of material statistics may be, and is 
often, overshadowed by the uncertainties inherent in design. 
Required statistics of structural steel are not generally 
available for common grades of structural steel because steel 
specifications and material specifications work with specified 
minimum values. Examining the existing literature on material 
properties of structural steel is, therefore, necessary and to 
obtain an estimate of the properties needed. Characteristic 
and representative sets of data were examined and estimates 
were made of the mean values and the coefficients of variation 
for tentative use. The principal material property affecting 
the resistance of a steel structure is the yield stress [22]. 
The values for use as proposed by T.V.Galambos et al is given 
in the table 3-3 overleaf. 



TABLE 3-3 Summary of Material Properties Used in LRFD 

Criteria 
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r — 

I Material Property 

Mean Value, in kips 
per square inch 

Coefficient of 
variation 

Modulus of 
elasticity in 
tension 

29000 

0.06 

Modulus of 
elasticity in 
compression 

29000 

0.06 

Modulus of 
elasticity in shear 

11200 

0.06 

Poisson's ratio 

0.30 

0.03 

Yield stress in 
flanges 

1.05 F y 

0.10 

Yield stress in webs 

1.10 Fy 

0.11 

Yield stress in 
shear 

0.64 F y 

0.10 


3.6 Variation of Safety Index 

The value of safety index may be varied to account for 
the importance of the structure. If the structure is 
important ( like public buildings, national monuments, places of 
worship, industries etc.), then they can be designed for a 
higher reliability factor, to take care of the stochasticity 
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of loading in these buildings. There are some structures in 
which failure of one or a few critical elements may result in 
the total loss of the structure ("weakest link" type 
structures ) in contrast to ductile or continuous structures 
("parallel" type structures) . 

Optimal levels of reliability for different types of 
structures could be obtained from an expected total cost 
optimization process. It could be decided, that 3=3.0 in 
ordinary buildings, 3=4.5 for very important buildings, and 
2 • 5 for temporary structures. It is possible to incorporate 
the statistical correlation between cross sections and between 
members and failure modes by suitably varying the value of the 
safety index 3 [10] . The LRFD formulation is versatile enough 
to incorporate these future developments in probabilistic 
design. 

3.7 Comments on LRFD codes 

The simple structures used by the LRFD approach to 
calibrate the load and resistance factors have closed form 
solutions, i.e., the response in these structures is available 
as an analytical expression in terms of the basic structural 
parameters. Therefore the limit state is also analytically 
available, making it easy to estimate reliability. However, 



45 


for most realistic structures the response is not available as 
a closed-form solution; it can only be evaluated through 
numerical procedures such as finite element analysis. 
Therefore, more complicated numerical procedures than those 
used in LRFD are needed to estimate the reliability of members 
in such structures. 

In the LRFD approach, however, the individual members in 
complicated structures are designed using the same load and 
resistance factors that were derived based upon the 
reliability analysis of simple structures. The use of isolated 
simple structures to derive safety factors is related to the 
basic design philosophy common to all codified design 
procedures. There are several advantages to the isolated- 
member approach: (1) In deterministic design methods that use 
factors of safety, preparing detailed requirements for each 
structural configuration is not practical; (2) the 
characteristics of the individual members and connections 
themselves are independent of the framework; and finally, (3) 
most research has been devoted to the study of such elements, 
and theoretical and experimental verification of their 
performance is readily available. Nevertheless, the 
performance of a member is directly dependent on its location 
in a structural configuration and its relationship or 
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connection with other members in the framework. Such 
dependence is not restricted to the computation of load 
effects through a deterministic analysis of the structure, but 
extends to the probabilistic characteristics of all the 
parameters of the structure. Only a probabilistic structural 
analysis of the entire structure can account for such 
influence and accordingly determine the risk or reliability of 
any individual member, thus enabling an improved approach to 
reliability-based design. 

An important objective of the reliability-based design 
methods such as LRFD is to reduce the scatter of nonuniform 
risk levels produced under various load combinations by the 
conventional design methods. As described in the AISC LRFD 
specification (Manual 1986) , the reliability indices inherent 
in the 1978 AISC specification (Manual 1978) , when evaluated 
under different load combinations and for various tributary 
areas of typical members, show a considerable range of 
variation. The LRFD approach seeks to narrow this range of 
variation of 3 values by specifying several "target" 3 values 
and selecting multiple load and resistance factors to meet 
these targets. Since the computation of 3 values in this 
approach is based on direct simulation using simple, isolated 
structural elements, an improved analysis of a realistic 
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structure would reveal that, even for the same load 
combination and the same limit state, there is considerable 
variation in the (3 values among the different members of the 
structure. Thus, there is further scope for improvement in the 
achievement of uniform risk by reducing the variation of (3 
values among different members in a structure, within the 
limitations of practical design. This is also advantageous 
from the point of view of structural optimization as in weight 
minimization, since uniform risk among members implies a 
balanced distribution of weight. 

A third aspect of the LRFD approach, which needs closer 
examination and possible improvement, is the consideration of 
the statistical correlation among the basic structural 
variables. The load and resistance factors in the LRFD 
approach were derived assuming statistical independence of 
variables. This may be reasonable for isolated simple members, 
which do not have too many variables, and assumption of lack 
of correlation may not significantly affect the determination 
of the reliability index. However, for members in structures 
such as frames, correlations among the random variables may 
have a significant effect, and so need investigation. 

The key to successful resolution of all these issues is 
the ability to perform reliability analysis of complicated 
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structures for which the response is not available as a 
closed-form solution in terms of the input variables, except 
in an algorithmic form such as finite-element code, like 
"NESSUS . " 



CHAPTER IV 


SYSTEM RELIABILITY ANALYSIS OF STEEL FRAMES 


4.1 Reliability Analysis of Complicated Structures 

Three types of solution strategies are possible for the 
reliability analysis of complicated structures; they are: 
(1) Direct simulation (2 ) approximation of the performance 

function by a polynomial; and (3) the stochastic finite 
element method [4]. 

The stochastic finite element method uses a more direct 
approach to the reliability analysis of structures. Starting 
with second-order statistics of the basic random variables, it 
keeps account of the variation of the quantities computed at 
every step of the deterministic analysis with respect to the 
basic random variables, and thus makes it possible to compute 
the statistics of response or the reliability for any limit 
state . 

For structures whose limit state is not available in 
closed form, Wu(1984) suggested the use of a simple, easily 
constructed second-degree polynomial that approximates the 
limit state in the neighborhood of the design point. Repeated 
deterministic analysis at selected points in the neighborhood 
and subsequent regression analysis are used to achieve this 
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objective. Then the Rackwitz-Fiessler algorithm is used to 
estimate the reliability index, through the solution of the 
approximate limit-state equation. 

Direct simulation, though robust is expensive. A large 
number of deterministic runs are required to compute the 
probability of failure, which is generally required to be very 
low in conventional civil engineering structures. The 
efficiency of the simulation can be improved by reducing the 
variance of the estimated probability of failure, which uses 
the same execution times and storage requirements without 
disturbing its expected value. Several such variance reduction 
techniques have been proposed and used in structural 
reliability analysis, e.g., importance sampling method. These 
variance-reduction techniques can also be combined further to 
increase the efficiency of the simulation. 

This chapter develops a method to quantify the effect of 
different types of collapse mechanisms of a structure under 
cumulative loading with the help of numerical examples. The 
purpose of the numerical examples is twofold: First, to 
illustrate reliability analysis of steel frames for the 
performance functions presented later in this chapter; and 
second, to examine steel frames designed using the LRFD 
approach and determine whether the target reliabilities of the 
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structures have been attained considering the overall 
structural configuration. 

4.2 Statistical Information 

For the reliability analysis, a probabilistic description 
of the variables is necessary. The stochastic variation of 
loads, material properties, sectional properties have been 
extensively studied in the earlier chapter; according to 
existing literature. Ellingwood et al. (1980) provided detailed 
statistical information, including the type of distribution, 
about some of these parameters. 

The dead load and all the resistance variables have been 
described as lognormal variables; wind load and the live load 
were described as type I extreme value variables. 

4.3 Statement of the Numerical Example 

Examine several steel structures, designed using the LRFD 
approach, without changing the structural configuration, but 
by varying the structural geometry and the loads. Determine 
whether the target system reliabilities, as stated by the 
codes are reached. Interpret the results obtained. 
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4.4 Analysis of the Problem 

One basic structural configuration for the plane steel 
frame is chosen. Using the same structural configuration, the 
structural geometry and the loads were varied to give sixteen 
different structures were obtained for analysis. The nominal 
values of the dead loads, live loads, wind loads are 
calculated to the best possible alternative, with the help of 
UBC codes [22] (1988). These values are tabulated for different 
runs, in table 4-1. 

4.4.1 Assumptions in analysis 

The following assumptions were made in the analysis: 

1. Elasto-plastic framed structures are used. If a moment 
exceeds the moment capacity at a section, a plastic hinge 
occurs and an artificial moment of magnitude equal to its 
resistant moment capacity is imposed at this section. 
Component failure due to buckling and violation of 
displacement constraints is not considered. 

2. The structural uncertainties are represented by 
considering only the moment capacities as random variables. 

3. Geometrical second-order and shear effects are 
neglected. The effect of axial forces on the reduction of 
moment capacities are also neglected. 
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4 . The order of loads and loading paths are not 
considered. 

These assumptions are often used in time-invariant system 
reliability analyses for ductile frame structures. 

4.5 Applied Plastic Design 

Until recent years most steel beams were designed based 
on elastic theory. The maximum load that a structure could 
support was assumed to equal the load that first caused a 
stress somewhere in the structure to equal the yield stress 
of the material. The members were designed so that the 
computed bending stresses for service loads did not exceed the 
design stress. Engineering structures have been designed for 
many decades by this method with satisfactory results. The 
design profession, however, has long been aware that ductile 
members do not fail until a great deal of yielding occurs 

the yield stress is first reached. This means that such 
members have greater margins of safety against collapse than 
the elastic theory seems to indicate. 

This sums up the basis of the plastic theory. The theory 
is that those parts of the structure stressed to the yield 
pdint cannot resist additional stresses. They instead will 
yield the amount required to permit the extra load or stresses 
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to be transferred to the other parts of the structure where 
the stresses are below the yield stress and thus in elastic 
range and able to resist increased stresses. Plasticity can be 
said to serve the purpose of equalizing stresses in cases of 
an overload. 

A statically determinate beam will fail if one plastic 
hinge develops. For a statically indeterminate structure to 
fail it is necessary for more than one plastic hinge to form. 
The number of plastic hinges required for failure of 
statically indeterminate structures will be shown to vary from 
structure to structure, but may never be less than two. 

One very satisfactory method used for plastic analysis is 
the virtual-work method. The structure in question is assumed 
to be loaded to its nominal capacity, Mn, and is then assumed 
to deflect through a small additional displacement after the 
ultimate load is reached. The work performed by external loads 
during this displacement is equated to the internal work 
absorbed by the hinges. For this discussion the small-angle 
theory is used. By this theory the sine of a small angle 
equals the tangent of that angle and equals the same angle 
expressed in radians. We can use these values interchangeably 
because the small displacements produce extremely small 
rotations or angles [ 19 ]. 
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This theory is the basis, of formulation of limit state 
functions in the analysis of the structures under 
consideration in the problem. 

4 . 6 Design of Structures 

A method is proposed to design the structures, in which, 
the objective is to estimate the element reliability under 
material nonlinearity represented by plastic hinge and not 
system reliability. The emphasis is on identifying the 
important linear segments of the nonlinear element reliability 
limit state through this procedure. In terms of 
implementation, the proposed method imposes a group of plastic 
hinges on the structure, instead of imposing only one hinge at 
each step as in current system reliability methods, as 
developed by Xiao, et al[7]. This grouped imposition is an 
important step that saves much computational effort for large 
structures since the number of structural reanalyses is 
greatly reduced. Particular group of plastic hinges, which 
will produce significant change, is isolated. 

4.6.1 Stepwise Design Procedure 

The algorithmic design procedure can be clearly seen in 
the flow chart as seen in Figure 4-1. 
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Figure 4-1 : Flow Chart for Design Procedure 
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The following steps were incorporated to design the structures 
under consideration: 

Step 1 - All the initial nominal values of the dead load D, 
live load L, and wind load W, are selected based on 
the UBC building codes, taking Nashville, TN, as the 
location center. The calculations* yielded the 
nominal values given in the table below. 

TABLE 4-1 NOMINAL VALUES OF LOADS 


Variable 

Nominal Value 

D - Dead Load 

4.0 kips 

L - Live Load 

7.5 kips 

W - Wind Load 

2 . 0 kips 


The basic structure to be analyzed is seen in Figure 
4-2, with dimensions and loading patterns. The live 
load and the dead load are applied at the center of 
each beam, with the wind load point application at 
the node of the column-beam junction. Based on this 
basic structural configuration, fifteen variations 
physically possible with variations in wind load, 




A 


10' (H) 


Y 


4 K (D)+7.5 K (L) 


4 K (D)*7.5 K (L) 
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Figure 4-2 : Basic Structural Configuration under analysis 


vertical dead and live load, horizontal bay 
dimensions and vertical height dimensions, were 
thought of; thus making a total of sixteen 
structures to be analyzed by the proposed method. 

The values of the wind load, live load, dead 
load, bay dimension, height dimensions are given in 
table 4-2 below for all 16 cases. 
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TABLE 4-2 VARIATIONS OF NOMINAL VALUES 


1 Case 

# 

Dead load 

Live load 

Wind load 

Bay size 

Height 



(D) ; kips 

(L) ; kips 

(W) ; kips 

(B) , feet 

(H) , feet 

Case 

1 

4.0 

7.5 

2.0 

18.0 

10.0 

Case 

2 

4.0 

7.5 

4.0 

18.0 

10.0 

Case 

3 

4.0 

7.5 

6.0 

18.0 

10.0 

Case 

4 

4.0 

7.5 

8.0 

18.0 

10.0 

Case 

5 

6.0 

10.0 

2.0 

18.0 

10.0 

Case 

6 

8.0 

12.5 

2.0 

18.0 

10.0 

Case 

7 

10.0 

15.0 

2.0 

18.0 

10.0 

Case 

8 

12.0 

17.5 

2.0 

18.0 

10.0 

Case 

9 

4.0 

7.5 

2.0 

18.0 

12.0 

Case 

10 

4.0 

7.5 

2.0 

18.0 

14.0 

Case 

11 

4.0 

7.5 

2.0 

18.0 

16.0 

Case 

12 

4.0 

7.5 

2.0 

18.0 

18.0 

Case 

13 

4.0 

7.5 

2.0 

20.0 

10.0 

Case 

14 

4.0 

7.5 

2.0 

22.0 

10.0 

| Case 

15 

4.0 

7.5 

2.0 

24.0 

10.0 

| Case 

16 

4.0 

7.5 

2.0 

26.0 

10.0 


Step 2 - The next step was to conduct the force study of all 
the 16 different structures. This was done with the 
aid of a structural finite element software, STAAD - 
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III, developed by Research Engineers Inc. The 
loading cases analyzed in this study according to 
LRFD[2] formulae were : 

U = 1.4D LRFD A4-1 (4-1) 

U = 1.2D + 1.6L LRFD A4-2 (4-2) 

U = 1.2D + 0.5L + 1.3W LRFD A4-4 (4-3) 

The outputs of the individual member forces were 
studied, with moment in the Z-direction being the 
prime governing factor as discussed in the 
assumption of analysis of the problem. The sections 
of maximum moment isolated, ford, for use in the next 
step of the design process that would be plastic 
design based on LRFD codes to design each member of 
all the 16 structures. 


Step 3 - A computer program* was developed on the lognormal 
distribution for 3 , as is discussed in Chapter III. 
The formula follows: 



(4-4) 


Where, 


R„= the mean resistance 



Q m = the mean load effects, 
which in our design process would be the plastic 
moment of the beam (including the effect of <|>)and the 
maximum moment induced in the beam derived from the 
force study. The target (3 and $ used for the columns 
and beams are given in the table 4-3 below: 


TABLE 4-3 TARGET VALUES OF (3 AND «J> FOR ELEMENTS 



3 

4> 

COLUMNS 

2.5 

0.85 

BEAMS 

__3^0__ 

0.90 1 


V r and V Q are the coefficients of variation of the 
plastic moment of the beam and the moment induced. 
Both of these coefficients of variations are 
numerically used as 0.10 as discussed in Chapter 
III. 

^fter using the program developed, all structures 
were designed for element reliabilities desired. 
The results of the design are tabulated in table 
4-4. 


* listed in APPENDIX 'A' 
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TABLE 4-4 DESIGN SECTIONS 




BEAM SECTION 

EXTERNAL 

INTERNAL 

CASE # 

(B) 

COLUMN SECTION 

COLUMN SECTION 



(Cl) 

(C2 ) 

1 

W8X28 

W8X10 

W6X9 

2 

W8X28 

W8X10 

W6X12 

3 

W8X28 

W8X15 

W8X15 

4 

W8X28 

W8X18 

W10X15 

5 

W10X33 

W10X12 

W6X9 

6 

W12X35 

W10X15 

W6X9 

1 7 

W10X49 

W6X25 

W6X9 

1 8 

W12X50 

W8X24 

W6X9 

| 9 

W10X22 

W6X12 

W6X9 

| 10 

W12X22 

W6X12 

W6X9 

11 

W12X22 

W6X12 

W4X13 

12 

W12X22 

W8X15 

W6X12 

13 

W8X31 

W6X15 

W6X9 

1 14 

W14X22 

W6X16 

W6X9 

I 15 

W10X30 

W8X15 

W6X9 

I 16 

W10X33 

W10X15 

W6X9 


The W-sections chosen for design are in accordance 
with the LRFD, AISC specifications. It should be 
noted that some sections are not practically 
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available in the market, but are listed in the 
codes. They are solely selected for theoretical 
purposes of close simulation to design conditions. 

Step 4 - A complete indeterminacy study was carried out of 
the structures. The following results were deduced: 
Number of possible hinges formed =10 

Number of redundancies = 6 


Number of independent mechanisms = 4 

Out of these 4 failure mechanisms, 3 critical 
mechanisms are identified- ( 1 ) beam mechanism; (2) 
column mechanism; (3) combined mechanism. 

Step 5 - The next step was the formulation of the 3 g- 
functions based on the virtual work study [23] of all 
the 3 mechanisms. The rotations, at the end of the 
plastic hinges are equal. The Figure 4-3 shows the 
different modes of failure possible in the 


structure . 



Mode 2 - Column Mechanism 



Note : All circlet indicate the formation and location of possile plastic hinges. 


Figure 4-3 : Significant Modes of Failure in structure. 
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The external work performed is always the 
product of the load and the average deflection of 
the mechanism. The average deflection equals one- 
half the deflection at the center or dominating 
plastic hinge. The internal work is the sum of Mp at 
each plastic hinge times the angle through which it 
works. The resulting expressions for the g- functions 
which is a result of difference between the internal 
resistance (work) and the external work performed. 
They are formulated below: 

1. BEAM MECHANISM: 

8.0*Mp b - (D+L) *B 

2. COLUMN MECHANISM: 

g 2 = ( 4 . 0*M pcl +2 . 0*Mp C2 ) - W*H 

3. COMBINED MECHANISM: 

g 3 =(2.0*M PC1 +3.0*M pC2 +12.0*M pb )-(W*H+B* (D+L) ) 

Step 6 - This step is explicitly explained in Chapter V, 
which concentrates on the formation of the fault 
tree risk analysis formation and discussion of the 


results . 



CHAPTER V 


SYSTEMS RISK ANALYSIS AND RESULTS 

System risk analysis is carried out using "NESSUS" by the 
development of fault trees that combine different modes of 
failure in the system. 

5.1 Fault Tree Analysis 

In calculating system reliability, it is important to 
include the probabilistic dependencies between multiple 
component failures, or between different failure modes. 
Failure to do so could result in significant errors. Fault 
tree analysis is a commonly used tool in risk assessment. A 
Fault tree is a mathematical construction of assumed component 
failure modes (bottom events) linked in series or parallel 
leading to a top event, which denotes the total system 
failure. A fault tree diagram essentially decomposes the main 
failure event (top event) into unions and intersections of 
subevents or combination of subevents. The decomposition 
continues until the probabilities of the subevents can be 
evaluated as single mode failure probabilities. The 
probabilistic fault-tree analysis is based on the limit state 
definition of the bottom events. Thus, one requirement for 
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system risk assessment is to compute failure function of each 
bottom event. Each bottom event is defined by a close form 
equation. 

A fault-tree has three major characteristics: bottom 

events, combination gates and the connectivity between the 
bottom events and the gates. The system risk assessment is 
limited to AND and OR gates. The OR gate implies that the 
output fault event is the union of subevents. The AND gate 
signifies that the output fault event is the intersection of 
the subevents. The different steps involved in the application 

of the fault-tree analysis method can be summarized as 
follows [24 ] . 

1. Development of a fault tree to represent the 

structural system. 

2. Construction of an approximate performance function 

for each bottom event. 

3. Determination of a dominant sampling sequence for all 

bottom events . 

4. Calculation of the system reliability using 

Adaptive Importance Sampling method. 

T° illustrate the Fault-tree analysis, consider a simple 
example consisting of two failure modes: yielding and 
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excessive displacement. Two failure functions can be expressed 
as, 

g, = R (Yield strength) - S (Stress) (5-1) 

g 2 = D (Allowable displacement) - d (displacement) (5-2) 

Failure occurs if [ g x <0 ] or ^9 <0] . Using standard 
probability notations, the system probability of failure is: 
Pf=P[(g,<0)u(g 2 <0)] (5-3) 

In general, 

Pf = P[g.<0]+P[g 2 <0]-P[( gl <0)n(g 2 <0)] = P,+P 2 -P 12 (5-4) 

In general, P 12 ranges from 0 to the smaller value 
of P x and P 2 therefore, P f ranges from [P 2 +P 2 ] to P 2 ( assuming 
P 2 >P X ) . Hence, by assuming independent events, the error 
ranges from -P 1 P 2 to P 1 (l-P 2 ). 

In application to the project, one OR gate is considered 
with three bottom events. The three bottom events represent 
the three failure modes of the structure. The representation 
of Fault-tree with three failure modes is shown in Figure 5-1. 
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SYSTEM FAILURE I 



x 

Combined Mechanism 


Figure 5-1 : Representation of fault tree 
analysis 


The Fault-tree analysis is carried out by two methods. They 
are: 

1. Adaptive importance sampling method. 

2. Standard Monte Carlo sampling method. 

5.1.1 Adaptive Importance sampling method 

Adaptive Importance Sampling is different from 
traditional importance sampling methods because of its ability 
to adjust automatically and by that reduce the sampling space. 
Because of this attribute, adaptive importance sampling method 
is highly efficient and accurate alternative for probabilistic 
analysis . 

Two options are available for selecting the sampling 
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boundaries. The first order adaptive sampling method uses 
hyperplanes, and the second-order adaptive sampling method 
uses parabolic surfaces. Both surfaces are constructed in the 
u-space and use the most probable point to define the 
beginning sample space. In general sampling space can be 
adjusted by increasing or decreasing the curvatures of the 
parabolic surface until there are no more failure points in 
the final sampling space. In the first order-based method, 
only the distance to the hyperplane is changed. In the second- 
order-based method, the curvature of the sampling boundary is 
updated first, then the final surface is shifted toward the 
origin [12] . 

5.1.2 Monte Carlo Sampling method 

Monte Carlo sampling method is a way of generating 
information for a simulation when events occur in a random 
way. It uses unrestricted random sampling (it selects items 
from a population so that each item in the population has an 
equal probability of being selected) in a computer simulation 
in which the results are run off repeatedly to develop 
statistically reliable answers. A sample from a Monte Carlo 


simulation is 

similar 

to a sample 

of 

experimental 

observations . 

Therefore, 

the results 

of 

Monte Carlo 
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simulations may be treated statistically. Monte Carlo methods 
are useful because they can handle very complex models, are 
guaranteed to work, and are exact in the limit as the number 
of samples becomes large. The disadvantage is that a very 
large number of simulations may be necessary [25] . 

5.2 Structural System Reliability using NESSUS 

System reliability considers failure at multiple 
locations, multiple failure modes, multiple components and 
combinations of all three. System reliability in NESSUS is 
currently addressed by a probabilistic fault tree 
analysis (PFTA) method. The driver module for system 
reliability is the SRA module with the PFTA methodology in the 
FPI module. The procedure implemented is intended to be 
accurate and efficient and build off the previous capabilities 
of NESSUS. 

The user defines system failure through the fault tree by 
defining the bottom events and their combination with "AND" 
and "OR" gates. Each bottom event considers a single failure, 
i.e., component reliability, and is defined through a finite 
element model and performance function. NESSUS will compute 
the reliability of each bottom event and a polynomial 
approximation, called a failure function, to the structural 
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response at the most probable point (MPP) using the AMV+ 
algorithm. The failure functions are then combined according 
to the fault tree [26]. System reliability is then computed 
using an adaptive importance sampling method. The adaptive 
importance sampling in this module has two features. First, 
the sampling region is focussed on the most important region 
where it has the highest probability of failure, and second, 
the sampling region is not predetermined. Instead, the 
sampling region is gradually increased by deforming the 
sampling boundary until the sampling region fully covers the 
failure region sufficiently. When the sampling region fully 
covers the failure region, the probability solution will 
converge, indicating that no more deformation is required [27] . 

There are several advantages to this approach. Because 
the failure functions are used, not just the probability of 
failure of each bottom event, the method can account for 
correlation between bottom events. The preexisting NESSUS 
capabilities for component reliability and failure function 
for each bottom event. In addition, adaptive importance 
sampling is typically an order or more faster than 
conventional Monte Carlo. 

The PFTA procedure implemented in NESSUS is being 
investigated for use with progressive fracture failure mode. 



5.3 RESULTS OF SYSTEM RELIABILITY CALCULATIONS 


The PFTA method uses the failure function about the MPP 
for each bottom event not just the probabilities. 

A summary of the system probalilities of failure and the 
respective safety indices for all the 16 cases are given in 
the table below 

TABLE 5-1 PROBABILITY AND SAFETY INDEX RESULTS 


1 Case # 

Probability of 

Safety Index ((3) 

1 

Failure (P f ) 


Case 1 

0 . 15728E-07 

5.534 

Case 2 

0.14167E-03 

3.630 

Case 3 

0 . 83901E-04 

3.763 

Case 4 

0.2191 6E-03 

3.516 

Case 5 

0. 56805E-10 

6.448 

Case 6 

0 . 54563E-12 

7.118 

Case 7 

0 . 27913E-12 

7.210 

Case 8 

0. 67462E-13 

7.401 

Case 9 

0 . 86963E-06 

4.782 

Case 10 

0. 85838E-05 

4.299 

Case 11 

0 . 47044E-04 

3.905 J 

Case 12 

0 . 51779E-06 

4.885 | 

Case 13 

0 . 97542E-09 

6.002 

Case 14 

0 . 21 091E-09 

6.246 
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Case 15 

0. 15191E-10 

6.645 

Case 16 

0 . 22340E-11 

6.922 


Considering the results, which we can even see in a 
graphical form as seen in figure 5-2, there is clearly a 
difference between the safety indices of the four different 
variations tried in the structural loadings and geometries, 
and we see that the system safety index is consistently higher 
than the target. This can be seen in the table 5-2 below; 


TABLE 5-2 SAFETY INDEX VARIATIONS 


CASE # 

AVERAGE SAFETY INDEX 

TARGET RELIABILITY 


( Pavg) 

INDEX (3 target) 

CASE 1 - CASE 4 

4.111 

3.0 

CASE 5 - CASE 8 

7.044 

3.0 

CASE 9 - CASE 12 

4.468 

3.0 

CASE 13 - CASE 16 

6.454 

3.0 
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These results are obtained by using the same LRFD 
criteria used for developing the codes, using first order 
probabilistic design. After these results were obtained, a 
sensitivity analysis was done for the first case, taking it as 
a representative structural configuration. The results of 
which are graphically represented in fig 5-3 to fig 5-5. From 
these graphs, we conclude that the live load is the most 
sensitive variable in the beam collapse mechanism and the 
combined collapse mechanism, and the plastic moment of the 
external columns is the most sensitive variable in the column 
collapse mechanism. 

It is decided to vary, the coefficients of variations of 
the vertical live load, plastic moment of the external column 
section, horizontal wind load, beam section; to study the 
effect of these variations of individual limit state variables 
on the safety index of the system. The results are depicted in 
fig 5-6 to fig 5-9. It is worth noting from fig 5-6, that even 
though the live load is the most sensitive parameter for two 
limit states, it has no effect at all on the system 
reliability. However, the plastic moment of the external 
columns, which is the most sensitive parameter for one limit 
state, affects the system reliability (safety index) 
considerably with change in its statistics. This means that 



77 


the plastic moment of the external column is the most critical 
parameter in the system. It also proves that even though a 
particular variable is sensitive in a single limit state, it 
may not be the most critical when system effects are 
accounted. 

5 . 4 Observations 

The main objective of this project - the validation of 
LRFD for actual structures - is achieved by comparing the 
reliability indices (computed using NESSUS) for the various 
limit states in a plane frame structure, with the target 
values used in LRFD. As seen in tables tabulate the NESSUS- 
computed p values for the structures and also the target 3 
values show that the former is consistently higher. 
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In cases 1-4, where, the design is determined by the 
variation of horizontal wind load, the p AVG is 4.111. In cases 
5-8 the design is dominated by the variations in vertical 
load, the P AVG is 7.044, which is more than twice the target p. 
The dominating factor in the design decision is variations in 
height of the structure, here the p AVG is 4.4 68. Finally the 
rest of the cases are dominated by change in bay dimensions of 
the structure, where the P AVG is 6.454. 

The NESSUS approach used here to estimate the system 
reliabilities designed according to LRFD differs from the 
latter in two respects: The effect of all structural variables 
is considered while estimating the reliability index for any 
P ar ticular case of the structural configuration whereas the 
LRFD method deals with the reliability of isolated members. In 
the approach described in this project, correlations are 
assumed between some random variables. Apparently, the load 
and resistance factors used by the LRFD approach are 
conservative, resulting in higher reliability of structures 
than the target reliabilities. The use of "standard" design 
situations such as simple beams, centrally loaded columns, 
tensile members, etc. to derive the load and resistance 
factors appears to have resulted in a conservative design for 
more complicated situations such as frames. 
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There is no appreciable effect of correlations between 
sectional properties on the reliability index. This is not 
surprising, since sectional properties show very little random 
variation (coefficient of variation is 0.1), whereas the 
largest variations are in the loading variables. Therefore, if 
there are any correlations among the load variables, it is 
possible that these might be more significant. 

If the proposed method, can give system structural 
reliability results to certain degree of accuracy, it is 
possible to use this method to formulate a procedure or 
relationship for optimum structural strength, ensuring uniform 
risk among different structural configurations. It is also, 
possible in the near future to relate the structural system 
reliabilities to the element reliability. 

However, it appears reasonable to account for the wind 
loading and enhance the value of yield stress. Also, it can be 
deduced that the two members are in two different 
configurations; therefore the combined effects of the random 
variables are different, altering the limit states and their 
distances from the origin. The observations also show that the 
reliability of a member is highly influenced by the structural 
configuration and that considering the effects of all the 
random variables is important. 



CHAPTER VI 


CONCLUSIONS AND FUTURE RESEARCH SUGGESTIONS 

6 . 1 Conclusions 

The study in this project covered the estimation of 
structural reliability under cumulative damage of single 
storey frame structures. The central idea was efficiently to 
impose the damage through a grouping operation, by exploiting 
the statistical correlations between modes of failures or by 
considering the amount of accumulated damage. Also it was 
sought to validate the load and resistance factors used in 
LRFD. The specific contribution of each finding is summarized 
and concluded as follows: 

The validation process involved comparison of the 
reliability levels achieved by the actual structure to the 
target reliability levels set up according to the LRFD 
cr ^ er ^- a * For numerical examples presented in this project, 
it is observed that the values of safety indices in the 
actual structure are higher than their target values. It is 
clearly seen from the fig 5-2, that if the framed structures 
designed according to the LRFD format, it leads to a 
higher safety index than is desired. This helps in concluding 
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that the structures designed with this method are 
conservative. Furthermore, the statistical correlations 
between the sectional properties of different members are 
observed to have an insignificant effect on the effect on the 
system reliability. The location of a member in a structural 
configuration appears to be much more significant. 

The LRFD method is based on the reliability analysis of 
simple, isolated members. Therefore, it does not consider the 
effect of the configuration on the stochastic response or 
reliability of a particular member; nor does it include the 
statistical correlations among all random variables of a 
structure. Using a finite element software like NESSUS offers 
the means to consider these factors and to estimate directly 
the reliability of the actual structural configuration. 
Therefore, this method can be used for a comprehensive 
validation of the LRFD approach, considering many different 
design situations. It is also possible to conduct sensitivity 
studies and compare the relative influence of various random 
parameters on the reliability. 

Finally the occurrence of non-uniform safety indices 
among different structural transformations, suggests that an 
approach with the assignment of sizes of critical members in 
each design group should be followed, based on the reliability 
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of the overall structure or safety index of the overall 
structure . 

6.2 Suggestions for Future Research 

Based on the results of the present study, the following 
topics may be addressed in future research: 

1. There is a need to use these results and similarly of 
many different structural configurations to determine the 
system affected element-reliability. In other words, the 
design factors may be derived for the system affected element 
reliability. 

2. The proposed study gives a very approximate estimation 
of system reliability under multiple time varying loads. The 
incorporation of the geometry (dimensions) of the structure as 
random variables would lead to a better result. 

3. The component resistances are assumed to be time- 
invariant. Practically the resistances are time varying due 
to aging, material deterioration etc. Fu and Moses evaluated 
the time dependent system reliability for a simple parallel 
system by updating the probability distribution of the element 
resistances at time t and using them to estimate system 
reliability at this time. However, the inclusion of time 
variant resistances in the estimation of system reliability 
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for a practical structure remains an important research topic . 
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.’HIS PROGRAM CALCULATES THE REQUIRED MEAN RESISTANCE , GIVEN THE 
1ETA VALUE AND THE LOADING EFFECT WITH THE COEFFICIENTS OF VARIATIONS 
)F THE RESISTANCE AND LOAD EFFECTS. 

-NITISH BERI 
26TH SEPT 95 

XCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

IMPLICIT REAL *8 (A-H,0-Z) 

WRITE (*,*) ' INPUT THE VALUE OF THE INDUCED MOMENT: (KIP-FT) ' 

READ ( * , * ) R 

WRITS (*,*) * INPUT THE VALUE OF C.O.V OF THE MOMENT:' 

READ ( * , * ) Vr 

' OPEN (UNIT=9, FILE= ' LNB . OUT ' , STATUS= ' UNKNOWN' ) 

) WRITE (*,*)' INPUT THE VALUE OF PLASTIC MOMENT OF SECTION-' 

READ ( * , * ) Q 

WRITE ( * , * ) ' INPUT THE VALUE OF C.O.V OF THE PLASTIC MOMENT-' 
READ ( * , * ) Vq 

OPEN (UNIT=9, FILE= ' LNB . OUT ' , STATUS= ' UNKNOWN' ) 

BETA1= (LOG (R/ Q) ) / ( (Vr**2+Vq**2) **0.5) 

WRITE (*, *) ' * 

WRITE (*,*)’ BETA = ', BETA1 

WRITE (*, *) ’ 

WRITE ( * / * ) ' IS THIS VALUE ACCEPTABLE ? ( 1=YES, 0=NO) ' 

READ ( * , * ) NUM1 
IF (NUM1 .EQ. 1 . 0) THEN 
GOTO 20 
ELSE 

GOTO 10 
ENDIF 

WRITE(9,*) ' + + *** + ** + ★****■*■*■* + + ■*•*** + ■*•**•*■****■*•**** + *** + + *.*. + ** f 

WRITE (9,*)’ RESULTS FOR THE BETA AND RESISTANCES OF MEMBER ' 
WRITE (9,*)' -NITISH BERI ' 

WR ITE (9, *) ' ********** + ****** + + *** + + *-*■* + ■*■ + ****•*■* + ** + **■■*•■*■** + + i 
WRITE (9, *) ' ' 

WRITE (9, 30) R 

WRITE (*, *) ' ********■*■ + + * + * + + ****•*•** + + ■*■•*■•*■* + + **** + + * + + + * + *•**** i 
WRITE (*, 30) R 

FORMAT ('THE VALUE OF MEAN RESISTANCE MOMENT= ',F5.2, ' KIP-FT') 

WRITE (*, *) ' ' 

WRITE (9,*)' 

WRITE (9, 40) Q 
WRITE (*, 40) Q 

FORMAT ( ' THE VALUE OF THE IMPOSED LOAD= ',F5.2,' KIP-FT') 

WRITE (*,*)’ ' 

WRITE (9, *) ' 

WRITE (9, 50) Vr 
WRITE (*, 50) Vr 

FORMAT ( ' THE VALUE OF C.O.V OF RESISTANCE = ' , F5.2) 

WRITE (*, +) » ' 

WRITE (9, *) ' ' 

WRITE (*, 60) Vq 
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WRITE (9, 60) Vq 

60 FORMAT (' THE VALUE OF C.O.V OF LOAD MOMENT^ ' ,F5.2) 

WRITE (*,*) ' 

WRITE (9,*)’ 

WRITE (*, 70) BETA1 
WRITE (9, 70) BETA1 

0 FORMAT (' THE VALUE OF BETA USED = ' , F5.2) 

T^TE ( * * ) 1 ** + **************'* , ****'* r ******* , '*'*'* r ******* , *****'*’ f 

STC 1' 

e:jd 
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loading calculations in accordance with the uniform building 

CODES : 


According to the Uniform Building Codes (UBC) ; 

1. Minimum roof live load = 20 psf. 

2. Minimum office live load = 50 psf. 

Intermediating these values, we use; 

Dead load on the frame = 20 psf. 

Live load on the frame = 40 psf. 

Choosing a bay width of 10 ft. 

Uniform dead load on beams = 20*10 = 200#/ft = 0.2 K/ft 

Uniform live load on beams = 40*10 = 400#/ft = 0.4 K/ft 

Therefore, 

Concentrated dead load at the center of beam = 3.6 K = 4 K 

Concentrated live load at the center of beam =7.2K=7.5K 

Basic Wind Speed = 70 mph. (Nashville) 

P=C e C q q s I 


1 = 1.0 


q s = 12.6 psf 
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C e = 1.06 (Exposure C) 

C q = 1.2 
Therefore, 

P= 1.06*1.2*12.6*1.0 
P= 16.03 psf. 

Linear load along column edge = 16.1*1 = 161#/ft = 0.161 

K/ft. 

Horizontal wind load = 1.61 K = 2 K 
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CASE 9 


4 K(D)+7.5 K(L) 


CASE 10 


4 K(D)+7.5 K(L) 


-J- 2 K(W)“ 


4 K(D)+7.5 K(L) 


CASE 12 


4 K(D)+7.5 K(L) 


y 


4 K(D)+7.5 K(L) 




4 K(DH7.5 K(L) 


CASE 14 


4 K(D)+7.J K(L) 


4 K(D)+7.S K(L) 



CASE 15 


4 K(D)+7.S K(L) 


4 K(D)-7.S K(L> 
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* Design Algorithm for running NESSUS for Probabilistic 
Analysis 

(Single limit state problems) 

1. Create a data file with a *.dat extension. 

2. Copy the *.dat file to for000.dat. This can be done by 
typing copy <*.dat> space for000.dat. 

3. To enter the failure function modify the subroutine 

respon . for 

4. To edit the file respon. for, type edt <respon.for>. 
This will create an editor asterisk on the screen. Type ' c' at 
the to get into full screen mode. 

5. Make changes and exit the file by holding the Ctrl key 
and pressing 'z' . This will again produce the editor's 
asterisk. Type 'save' at the asterisk and close the file. 

6. Once the subroutine is modified, it has to be compiled 
and linked to the library. To compile type fortran 
<f ilename . f or> 

7. Link the compiled file to the library by typing link 
filename (omit extension), nes/lib 

8. The probabilistic analysis can be done by typing run 


nessus at the VAX prompt. 
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9. NESSUS will then ask for the filename. The filename 
should be typed without the .dat extension. 

10. Once the run is completed, the output information 
will all be stored in for000.dat file. To preserve this 
information, change the name of the for000.dat to <input 
f ilename . out>, by typing ren for000.dat space <input 
filename . out> 

11. To study the sensitivity analysis results, type 
<input filename.mov>. if the safety index is low or the 
probability of failure is high, identify the most sensitive 
design parameter from the sensitivity analysis. 

12. Increase/decrease the coefficient of variance of the 
most important design parameter and do the probabilistic 

analysis again. This can be done by repeating steps 6 through 

11 . 

* Design Algorithm for Running NESSUS for Probabilistic 
Analysis ( Multiple Limit state problem / System reliability) 

1. Create a data file with a *.dat extension. 

2. Copy the *.dat file to for000.dat. This can be done by 
typing copy <*.dat> space for000.dat. 

3. The probabilistic analysis can be done by typing run 
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nessus at the VAX prompt. 

4. NESSUS will then ask for the filename. The filename 
should be typed without the .dat extension. 

5. Once the run is completed, the output information will 
all be stored in for000.dat file. To preserve this 
information, change the name of the for000.dat to cinput 
filename . out>, by typing ren for000.dat space Cinput 
filename . out> 

6. Increase/decrease the coefficient of variance of the 
most the most sensitive design parameter and do the 
probabilistic analysis again. This can be done by repeating 
steps 1 through 5 . 
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CHAPTER I 


INTRODUCTION 


1-1. Background 

In modern engineering design, the need for designing 
high-reliability, optimal structure systems has been 
increasing recently due to the demand for greater quality, 
reliability, and lower cost or weight. An optimal structure 
system design must satisfy the performance of cost, volume, 
weight, or speed ratio objectives as well as the system or 
component reliability constraint. The latter is used to 
quantify the uncertainty existing in different failure models, 
loading conditions, material properties, and geometric 
parameters. To deal with these uncertainties, reliability 
technology provides tools for formal assessment and analysis. 
Meanwhile, optimization technology plays an important part 
to meet the optimal design objectives. Therefore, the 
combination of reliability and optimization technologies is 
a viable way to design high-reliability optimal structural 
systems . 

The scope of this study is, however, mainly concentrated 
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on the design of a gear train using reliability-based 
optimization design method. Some failure models of gear, such 
as wear and thermal conditions, are not investigated in this 
study. 

1-2. Objectives and Organization of the Report 

2.1. The theory of the probabilistic design methodology in 
depth and an overview of the software, NESSUS (Numerical 
Evaluation of Stochastic Structure Under Stress), are 
described in Chapter II. 

2.2. The theory of GRG (Generalized Reduced Gradient method 
for optimization) and its application, and development of GRG 
computer program for this project, are described in detail in 
Chapter III. 

2.3. Discussion of the combination of reliability design 
method and optimization design method is presented in Chapter 
IV. 

2.4. The application of reliability-based optimization design 
method for a gear train and the results, their interpretation, 
explanation and comparison are described in Chapter V. 

2.5. The summary and conclusions of the present study 
and suggestions for future research are presented in Chapter 


VI. 



CHAPTER II 


PROBABILISTIC DESIGN METHODOLOGY 

2-1 Introduction 

In engineering designs, decisions are often required 
irrespective of the state of completeness and quality of 
information, and thus are made under conditions of 

uncertainty. In other words, the consequences of a given 
decision cannot be found out with complete confidence. Besides 
the fact that the information must often be inferred from 
similar circumstances or derived through modeling, many 
problems in engineering involve natural processes and 
phenomena that are inherently random. The states of such 
phenomena are naturally indeterminate and thus cannot be 
described definitely. For these reasons, decisions required 
during engineering design invariably must be made under 
conditions of uncertainty. 

The effects of such uncertainty in design are important. 
To be sure, the quantification of such uncertainty and 
evaluation of its effects on the performance and design of an 
engineering system should include concepts and methods of 
probability. Furthermore, under conditions of uncertainty, 
the designs of engineering systems involve risks, and the 
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formulation of related decisions requires them to be risk- 
free. The problems of uncertainty in the design can be 
overcome by applying the methods of probabilistic design. 
Thus, the role of probability is quite pervasive in 
engineering. It ranges from the description of information to 
the development of bases for design and decision-making [1]. 

PDM (Probabilistic Design Method) is concerned with the 
probability of non-failure performance of structures or 
machine elements. It is much more useful in situations in 
which design is characterized by complex geometry, possibility 
of catastrophic failure, or sensitive loads and material 
properties. Current studies on probabilistic structure 
analysis methods have resulted in a new class of tools that 
the engineer can use to obtain direct information on the 
uncertainty of structural performance. Reliability analysis 
evaluates the probability by a rational treatment of the 
uncertainties in various design parameters. It is becoming 
substantially evident that the PDM is beginning to attract 
more attention [2]. The PDM has been successfully applied to 
various loading conditions encountered during space flight 
[3]. Some reasons for the increasing acceptance of the PDM [2] 
are 


1) The deterministic method can provide some basic 
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information to complex design problems but provides no 
information with regard to the reliability of the design. 

2) Probabilistic computations are becoming simpler and less 
expensive because of software being developed. 

3) The PDM and the information it provides are becoming more 
widely understood and better appreciated. 

One of the most recent computer codes is NESSUS 
(Numerical Evaluation of Stochastic Structure Under Stress). 
This code was developed under NASA's probabilistic structure 
analysis program. An overview of NESSUS and the description of 
its development are given by Cruse et al.[4,5]. 


2-2 . Application of PDM 

Because probabilistic design method (PDM) is concerned 
with the probability of failure-or preferably, reliability-it 
is most useful when uncertainties in material properties and 
loading conditions are considered. To apply probabilistic 
design methodologies (PDM) , all uncertainties are modeled as 
random variables, with selected distribution types, means, and 
standard deviations. The primitive (random) variables that 
affect the structural behavior have to be identified. Every 
design project demands some sequential stages of reflection 



6 


before one can arrive at the final design goal. This is also 
the case with PDM. The various design stages of PDM are as 
follows : 

1) Defining the problem. 

2) Generating design parameters. 

3) Relating the defined problem to the design 
parameters . 

4) Assembling data and applying probability 
concepts . 

5) Using probabilistic Analysis. 

6) Interpreting results. 

The design stages of PDM are shown in Figure 2-1. 



7 



Figure 2-1: Design stages in PDM [6] 
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2-2.1. Problem definition 

The first step which a designer takes in solving a design 
problem is to find out the main objective of the design. After 
finding out the objective, the next step is to define in a 
precise manner the functional requirements of the system or 
component to be designed. These functional requirements should 
be able to characterize completely the design objective by 
defining it in terms of specific needs. With a clear 
understanding of what one is searching for, the designer then 
goes to the next stage. 

2-2.2. Generating design parameters 

In order to solve the defined problem, acceptable design 
parameters must be generated that will meet the defined 
functional requirements. To generate the design parameters, 
one uses an appropriate design model. The various parameters 
(loads, material properties, geometry, etc.) are taken into 
consideration. The design parameters to be selected depend on 
the objective of the design [6] . 

2-2.3. Relating the defined problem to the design parameters 

After defining the design parameters, the designer then 
relates the functional requirements in the functional domain 
to the design parameters in the physical domain, to be sure 
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that the objective is satisfied. If the relation is 
satisfactory, the designer goes to the next stage. If the 
relation is not satisfactory, it is redefined, so that the 
objective is satisfied. 

2-2.4. Data assembling and application of probability concepts 

This stage requires assembling the essential data that 
are available on the problem with regard to the design 
parameters. If some data are unavailable, then it becomes 
necessary to perform a computational simulation analysis to 
generate the missing details. Once the data have been 
assembled, the next stage is to analyze the assembled data. 
NESSUS is the computer tool used to perform the analysis. 
NESSUS has three modules, known as NESSUS/PRE, NESSUS/FEM, and 
NESSUS/FPI. 

NESSUS/PRE is a preprocessor, which prepares the 
statistical data needed for the probabilistic design analysis. 
It allows the user to describe the uncertainties in the 
structural design parameters. The uncertainties in these 
parameters are specified by defining the mean value, the 
standard deviation, and the distribution type, together with 
an appropriate form of correlation. Correlated random 
variables are then decomposed into a set of uncorrelated 
vectors by a model analysis. 
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NESSUS/FEM is a general purpose finite element code, 
which is used to perform structural analysis and evaluation 
of sensitivity due to variation in different uncorrelated 
random variables. The response surface, defined in terms of 
random variables required for probabilistic analysis in 
NESSUS/FPI, is obtained from NESSUS/PRE. NESSUS/FEM 
incorporates an efficient perturbation algorithm to compute 
the sensitivity of random variables [6] . 

NESSUS/FPI is an advanced reliability module, which 
extracts the database generated by NESSUS/FEM to develop a 
response model in terms of random variables. In this module, 
the probabilistic structural response is calculated from the 
performance model. The probability of exceeding a given 
response value is estimated by a reliability method. Inside 
the NESSUS/FPI module is a sensitivity analysis program, which 
determines the most critical design parameters in the design. 
The input data for NESSUS/PRE requires fundamental knowledge 
of statistics or probability theorems. The expected details 
will include determining the mean, standard deviation, median, 
coefficient of variation, variances, etc., associated with 
each random variable. The designer also determines the 
probability distribution function that best describes each 
random variable. The different modules of NESSUS are shown in 


Figure 2-2. 
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2-2.5. Probabilistic Analysis 

It is at this stage of the design that the designer 
defines a limit state function. The limit state function is a 
function that defines the boundary between the safe and 
failure regions. In the limit state function approach for 
structural reliability analysis, a limit state function g (X) 
is first defined. The g-function is a function of a vector of 
basic random variables, X =(X x , X 2/ X 3 , . . .X n ) with g(X)=0 

being the limit state surface that separates the design space 
into two regions, which are the failure g(X)<0 region and 
the safe g (X) >0 region. Geometrically, the limit state 
equation, g(X)=0, is an n-dimensional surface that may be 
called the "failure surface." One side of the failure surface 
is the safe state, g(X)>0, whereas the other side of the 
failure surface is the failure state, g(X)<0. 

The probability of failure in the failure domain Q is 
given by 


P f = Jo* • *Jf*(X) dx (2-1) 

where f x (X) is the joint probability density function of X and 
Q is the failure region. The solution of this multiple 
integral is, in general, extremely complicated. Alternatively, 
a Monte Carlo solution provides a convenient but usually time- 
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consuming approximation. The limit state function method uses 
the Most Probable Point (MPP) search approach, shown in Figure 
2-3. The Most Probable Point is the key approximation point 
for the FPI analysis; therefore, the identification of MPP is 
an important task. In general, the identification of the MPP 
can be formulated as a standard optimization problem and 
solved by proper optimization methods. 

From the Figure 2-3, as the limit state surface, g(X)=0, 
moves closer to the origin, the safe region, g(X)>0, decreases 
accordingly. Therefore, the position of the failure surface 
relative to the origin of the reduced varieties should 
determine the safety or reliability of the system. The 
position of the failure surface may be represented by the 
minimum distance from the surface g(X)=0 to the origin. The 
point on the surface with minimum distance to the origin is 
the Most Probable Point (MPP) . This is usually determined by 
fitting a local tangent to g (X) and moving this tangent until 
MPP is estimated. 

In the NESSUS code, MPP is defined in a transformed space 
called u-space where the u's are independent to simplify the 
probability computations. By transforming g(X) to g(u), the 
most probable point, u*, on the limit state, g(X)=0, is the 
point that defines the minimum distance from the origin to the 
limit state surface. This point is most probable (in the 
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u-space) because it has maximum joint probability density on 
the limit state surface. The required minimum distance is 
determined as follows. The distance from a point u*=(u 1 * / u 2 ‘, 
..., u n *) on the failure surface g(u)=0 to the origin is 

" i - 

D - ( £ “i > 2 (2-2) 

i - 1 

where D is the minimum distance from the point on the limit 
state surface to the origin. 

The FPI code assumes only one MPP. In general, however, 
the possibility exists that there may be multiple local and 
global Most Probable Points. A two MPP problem can occur; 
for example, if the g-function is quadratic, the search 
algorithm may result in oscillating (non-convergent ) search. 

Several approaches are available to search for the MPP. 
The search procedure depends on the forms and the number of 
the g-function (s) . One efficient method in use is the Advanced 
Mean Value method (AMV) . This method blends the conventional 
mean value method with the advanced structural reliability 
analysis method. This method provides efficient cumulative 
density function analysis and the reliability analysis. The 
step-wise AMV method can be summarized as follows [7] : 


1. Obtain the g (X) function based on perturbations about 
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the mean values. 

2. Compute the cumulative density function of the 
performance function at selected points using the 
fast probability integration method. 

3. Select a number of cumulative density function values 
that cover a sufficiently wide probability range. 

4. For each cumulative density function value, identify 
the most probable point. 

The analytical process involved in the limit state 
approach can be illustrated by a basic structural reliability 
problem. In the problem, only one load effect, S, limited by 
one resistance, R, is considered. 

If one considers a case when R and S are independent, 
the limit state equation can be expressed as 

g = R - S (2-3) 

and the probability of failure can be expressed as; 


P r P(R - S s 0) • /“/“/*( r )// s) dr ds 


(2-4) 
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For any random variable the cumulative density function F(x) 
is given by 

F x (x) = P(Xsx) = J4(Y ) dy (2-5) 

provided that x z y. Therefore P f is expressed as 

P f = P(R-S^O) = / F R ( x ) f s ( x ) dx (2-6) 

Assuming a special case of normal random variables, for 
some distributions of R and S, it is possible to integrate the 
equation (2-6) analytically and find out the probability of 
failure. If S and R have mean y R and p s and variance a R 2 and a s 2 
respectively, the g-function has a mean p g and variance a g 2 , 
given by 


P g = Mr - Ms 

(2-7) 

o g 2 - ° R 2 + o s 2 

(2-8) 


Therefore, the probability of failure is given as 


P f - P( R - S s 0 ) = P( g <: 0 ) . 4>[ 


0 - 


»*r 


o 

t 


] 


(2-9) 
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Which reduces to: 


* [ - - p == = l ] ■ * (-P) (2-10) 

V< ) 

where 3 is defined as the safety index. 


P ' — (2-11) 

a 

t 

Thus, the probability of failure is given as 

P f - 4> (-P) (2-12) 

which can be written as 


P f - 1 - 4>(P) ( 2 - 13 ) 

Reliability is the probability that the structure will not 
violate a given performance criterion during a specified 
period. This can be mathematically expressed as 


Pr - 1 Pf 


( 2 - 14 ) 


where P r is the reliability and P f is the probability of 
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failure. Structural reliability analysis evaluates the 
probability of failure by rationally treating the various 
uncertainties . 


2-2.6. Interpretation of Results 

This is the last stage in the methodology. When the 
designer approaches this stage, one interprets the results 
obtained about the initial objective. If the results do not 
satisfy the functional requirements in the stage 1, the 
designer may adjust design parameters to achieve the set 
objective . 


2-3. Probability Sensitivity Factors 

In engineering performance analysis many sensitivity 
measures can be defined. Knowing the effect of each random 
variable in the analysis is important for the designer. The 
sensitivity information is quantified by sensitivity factors. 
Sensitivity factors suggest which random variables are crucial 
and require special attention. 

The commonly used sensitivity factor in deterministic 
analysis is the performance sensitivity, dZ/dX. ± , which measures 
the change in the performance due to the change in a design 
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parameter. This concept can be extended to the probabilistic 
analysis in which a more direct sensitivity measure is the 
reliability sensitivity that measures the change in the 
probability/reliability relative to the distribution 
parameters such as the mean and the standard deviation. 
Although not automated in the code, this analysis can be 
performed by varying the parameters. 

Another, perhaps more important, kind of probability or 
reliability sensitivity analysis is the determination of the 
relative importance of the random variables. This analysis can 
be done, for example, by repeated probabilistic analysis in 
which one random variable at a time is treated as a 
deterministic variable. The results of the analyses, for 
example, are a number of cumulative density function curves or 
reliabilities. Based on the results, the relative importance 
of the random variables can be analyzed. The standard FPI 
output includes a first-order sensitivity factor that provides 
approximate relative importance of the random variables. 



CHAPTER III 


OPTIMIZATION DESIGN METHODOLOGY 


3-1 . Introduction 

Optimization is the method of obtaining the best result 
under given circumstances. In design, construction, and 
maintenance of any engineering system, engineers have to take 
many technological and managerial decisions at several stages. 
The ultimate goal of all such decisions is either to minimize 
an effort required or to maximize a desired benefit. 
Engineering design is a multiphase process requiring constant 
decision making by the designer. Based on his decision, the 
engineer is able to define variables, a design objective, and 
a set of constraints that must be met in order that the design 
is a workable solution. By developing corresponding equations, 
the design problem can be formulated into a standard form 
acceptable to mathematical programming techniques. This 
standard form is defined below. 

Minimize f(X) 

X = (x 1 ,x 2 ,x 3 , . . .,x N , ) T , X e R N (3-1) 
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subject to 


4>k(X ) * 0 k = 1,2,3,..., K (3-2) 

*,(X) = 0 1 = 1,2,3,...,L (3-3) 


where 

X is a column vector of design variables 
N is total number of design variables 
f (X) is objective function 

0 k (X)is K inequality constraint functions 
(X) is L equality constraint functions 


A more general occurrence in engineering design arises 
when expressions (3-1 to 3-3) are nonlinear. This is known as 
the nonlinear programming (NLP) problem. No general method has 
been developed to solve nonlinear problems in the sense that 
the simplex algorithm exists to solve the linear problem. 
Although many strategies have been suggested, comparative 
studies [8,9] have shown that no method has been successfully 
applied to all problems. In this project the Generalized 
Reduced Gradient (GRG) method will be described. The GRG 
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method avoids many of the problems associated with penalty 
function and LP-like methods [10], producing one of the most 
powerful methods currently known for handling the constrained 
nonlinear programming problem. The principle behind this 
method is quite simple, but its application is rather complex 
[ 11 ] . 


3-2. Theory of Generalized Reduced Gradient Method 

The Reduced Gradient method was originally given by Wolfe 
for a nonlinear objective function with linear constraints 
[12,13]. A generalization of Wolfe's method to accommodate 
nonlinearities in both the objective function and constrains 
was first accomplished by Abadie [14]. Concurrently to both 
Wolfe and Abadie, Wilde and Beightler developed their 
differential algorithm based on the constrained derivative 
[15] . The constrained derivative and the reduced gradient 
employ much the same theoretical basis, but for purposes of 
this discussion, the method shall be known as the Reduced 
Gradient method. The case of nonlinear constraints was 
pioneered by Abadie [14], who called it the "generalized 
reduced gradients (GRG)." Later variants were developed by 
Lasdon and Waren [16], Gabriele and Ragsdell [10]. Both 
Gabriele [17] and Lasdon and Waren have implemented versions 
for large sparse systems. The general constrained nonlinear 
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programming can be stated in the following form: 


Minimize f(X) 

X = (x„ x 2 , x 3) ..., x N , ) T , x e R N (3-4) 


subject to 


f m (X)= 0 m = 1, 2, .... M 


(3-5) 


AsXiB 


(3-6) 


The N x l vectors A and B represent upper and lower bounds on 
the design vector X. These upper and lower bounds can be 
assumed to be the finitive or infinitive bounds. The 
inequality constraints have been included as equality 
constraints by using the following transformation: 


\MX) = <M X ) - S k = 0 

0 * S k s ~ k = 1, 2, . . . / K (3-7) 
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If 4> ( X )* 0, the equation (3-7) will be changed to 
iM X ) = <M X ) + S k = 0 

— £ S k <; 0 k = 1, 2, . . . , K (3-8) 

The variables S k are slack variables that are included in the 
original set of design variables. Therefore, the parameter N 
represents the total number of design variables plus the 
number of slack variables used for the transformation of 
(3-7) or (3-8) . The parameter M represents the total number of 
constraints : 


M = L+ K 

Where 

L is number of equality constraints 
K is number of inequality constraints 


(3-9) 


It should be stressed that the constraints of (3-5) are 
nontrivial constraints; that is, they are functional 
constraints. Variable bounds are defined in (3-6) and will 
require separate handling. 

Linearization of equation (3-5) will result in M equality 
constraints with N independent variables. If the constraints 
were linear, all we have to do is to use elimination process 
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to reduce the number of independent variables to K using the 
equality constraints and then substituting the independent 
variables into the objective function f ( x ). Unfortunately, 
the problem is nonlinear so direct substitution is very 
difficult. Consider the following strategy whose fundamentals 
can be found in the simplex method of linear programming. 
Divide the design vector of equation (3-9) into two classes 
that shall be known as the decision and state variables. 


X = [ 

z , Y ] T 

(3-10) 

z = [ 

^1 / ^2 / • • • ] 

(3-11) 

Y = [ 

Yl/Y; 2/ • • - Ym ] t 

(3-12) 


where 

Z : decision variables; y : state variables. 

Q : number of decision variables, Q = N - M. 

The decision variables are completely independent, and the 
state variables are slaves to the decision variables used to 
satisfy the constraints \|/ ( X ) . 

The following notation will be useful in the discussion 
to follow: 


g(T) 


wo t WO 

dy x ’ dy 2 


WO ,r 


(3-13) 
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g(Z) - 


dm ' dm 

dz, ’ dz 2 


dm 

dz„ 


V 


( 3 - 14 ) 


dfr, a^, 

dz x dz 2 


d± 

dZ 


d *u d * 


u 


dz, dz , 



dijr 


u 


dz 


Cj 


<W «0 


( 3 - 15 ) 


a *, ^^1 

dy 2 

djr 

ar " 

ay, ay 2 


aj, 

ay w 


a+ 


A/ 






( 3 - 16 ) 


Let us examine the first variation of f (X) and (X) , 


df - g(Z) T dZ ♦ g(V) T dY 


( 3 - 17 ) 
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dty - — dZ + ^ dY * 0 
BZ BY 


(3-18) 


where 

dZ = Q* 1 vector of differential displacements of Z 
dY = M*l vector of differential displacements of Y 

Solving (3-17) and (3-18) and rearranging will yield the 
following linear approximation to the reduced gradient: 


dY . .**li± d z 

BY BZ 


(3-19) 


Substituting (3-19) into (3-17) 


g r (X) T - g(Z) T - g(Y) 


t d ^' 1 

BY 


BZ 


(3-20) 


The reduced gradient defines the rate of change of the 
objective function with respect to the decision variables with 
the state variables adjusted to maintain feasibility. 
Expression (3-19) gives the changes necessary in the states 
for a given change in the decisions for linear constraints. 
Geometrically the reduced gradient can be described as a 
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projection of the original N - dimensional gradient onto the 
( N - M ) - dimensional feasible region described by the 
decision variables. 

A necessary condition for the existence of a minimum of 
an unconstrained nonlinear function is that the elements of 
the gradient vanish. Similarly, a minimum of the constrained 
nonlinear function occurs when the appropriate elements of the 
reduced gradient vanish. This conclusion can be verified by 
a comparison with the Kuhn - Tucker [18] conditions for the 
existence of a constrained relative minimum. 

By first transforming the variable bounds into inequality 
constraints, 


$ i (X)= x L ~ ai *0 (3-21) 

<t> i+N (X)= bj. - Xi *0 (3-22) 


where i = 1,2,3, ...,N 

we can form the following lagrangian function: 


L(X,U,W) - f[X) 


U 2 N 

E w . ♦.(*) - E u,*/*) 

m-l >1 


( 3 - 23 ) 


The following Kuhn-Tucker necessary conditions hold for a 
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point to be a relative minimum X*, 


g(X‘) T ♦ E 

K 

dilr J di>. 

" -T u; ' - o 

(3-24) 

M-l 


dX j. i d(X) 


and 




Ur m (X*)= 0 


m = 1, 2, . . . ,M 

(3-25) 

$5 (X* ) * 0 


j = 1, 2, . . . , J =2N 

(3-26) 

Uj* <t>j(X*) = 

0 

j = 1, 2, . . . , J =2N 

(3-27) 

U a * i 0 


j = 1/2/. . . f J = 2N 

(3-28) 

w m * * o 


m = 1/ 2, . . . ,M 

(3-29) 

Introducing decision 

and state variables into 

(3-24) 

decomposing we can obtain the following form: 


g ( Z‘ ) r * W' T 

dijr 

■ - U' T d ^ =0 

(3-30) 


dZ' 

az- 


g ( v ) T . w T 

di| r 

. - w T 8,|> . 0 

(3-31) 


dY‘ 

dY‘ 



For the reasons that a state variable is not allowed to 
be equal or sufficiently close to either of its bounds. Form 
expression (3-27) the elements of U* corresponding to the 
state variable bounds must be zero. Also, those elements of 



31 


3<t> /dY corresponding to the decision variable bounds will be 
zero eliminating the last term of (3-31) . Solving (3-31) for 
W* and substituting into (3-30) will produce the following 
expression: 


* ( z- ) r • g ( r- ) r it 

ar- 


i*. . D- T iS- . 0 

dZ' dZ‘ 


(3-32) 


Rearranging (3-32), we obtain 


i* . , ( Z- / - * ( F- )' it! i* (3-33) 

dZ' BY' dZ' 

It can be recognized that the right-hand side of (3-33) 
is g r (X) . By examining the possible values of the left-hand 
side of (3-33), a candidate point X will be X* if 

g r (X)i > 0 if Zj = a t 
g r (X) t < 0 if z A = bi 
g r (X) i = 0 if s Zj s bj 


where 


i 


(3-34) 
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3-3. Generation of Feasible Starting Points 

In the application of GRG method, generation of feasible 
starting points is an important step. As the final element of 
problem presolution analysis, the formulation should be tested 
for feasibility. While the preceding stages of problem 
preparation and analysis may have resulted in a numerically 
stable, bounded, and nonredundant formulation, it is always 
possible that, as a result of poor data or calculation errors, 
the problem constraints simply exclude all possible solutions. 
Thus, whether or not the optimization algorithm selected for 
use requires a feasible starting point, it is good practice to 
devote some effort to generating a feasible starting point. 
Obviously, if no feasible starting point can be obtained, 
there is little point in proceeding with optimization. 
Instead, the model must be inspected again and validated in a 
piecemeal fashion until the sources of error are identified. 
If a feasible point can be generated, and if the variable 
values at the generated point appear reasonable, then one can 
proceed with problem solution with a fair degree of confidence 
that the optimization runs will be productive. 

A very common way of generating feasible starting points 
is direct minimization of the constraint infeasibilities. The 
method of minimization of unconstrained penalty-type functions 
has been proposed. This method is often preferable. 
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especially, for higher dimensionality and tightly constrained 
problems [19] . The procedure consists of solving an 
unconstrained minimization problem whose objective function is 
an exterior penalty function. The starting point is thus 
obtained as the solution of the problem. 

Minimize : 


k 


N 


A * 0 * E (*, (*)) a ♦ E ( Q A C *))) 2 

' ■ 1 i - 1 


(3-35) 


where 

(X) = quality constrained functions, 
k = number of quality constrained functions. 

(X) = inequality constraints of the variable bounds. 

N = number of unknown variables in function. 

Clearly a feasible point is one that will result in an 
objective function value of zero. Hence, the unconstrained 
minimization will be terminated when f (X) becomes sufficiently 
small. Generally, the minimization can be simplified if the 
problem is posed in equality-constraint-free form. 

3—4 . Perform the Line Search to Locate Local Minimum 

When a search direction D (Y) is determined for the state 


34 


variables, the vector D(Z) has defined a line in the reduced 
N - M dimensional space along which exists a local minimum of 
the objective function. It is the task of this section to 
locate the minimum so that it might be used as a starting 
point for the next iteration. The performance of this task is 
a common occurrence in many unconstrained searching 
techniques, and in the case of the reduced gradient represents 
the bulk of the computational effort. 

The normal course of events in locating a minimum along 
a line consists of two phases. The first phase involves 
locating an initial bracket within which the minimum is known 
to be contained. This is commonly referred to as the bounding 
phase. The second phase would consist of some efficient scheme 
of narrowing the initial bracket unit the minimum is known to 
be within some tolerance. Both these phases are outlined in 
more detail in [20,21,22]. 

The Reduced Gradient Method uses this same two - phase 
procedure with modifications to accommodate the use of state 
and decision variables. From a starting point (Z,Y) k we move to 
a new point (Z,Y) k+1 according to the step prescription 

Zi k+1 = b ± if Zi k + a D(z ± ) * b ± 

Zi k+1 = ai if Zi k + a D(z ± ) s a A 
z ± k+1 = z L k + ck D(Zi) otherwise 
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where i = 1,2,3, ...,Q (3-36) 
and 

Y m k+1 = Y m k + a D(YJ m = 1,2, . . . ,M (3.37) 
where a = step length parameter. 

Because of nonlinearities arising in the constraint 
functions, the point (Z,Y) k+1 is likely to be infeasible. 
Holding the decision variables Z k+1 constant, the state 
variables Y k+1 are adjusted to obtain a feasible point, 
(Z, Y) k+1 . This situation is shown in Figure 3-1 for the case 
M = 1/ Q = 1. This step is equivalent to the solution of M 
nonlinear equations (i|/(X)* 0) in M unknowns (Y) . A number of 
numerical techniques are available in the literature to 
perform this task. Newton's method [23] has proven to be an 
efficient technique as well as convenient since the necessary 
partial derivatives have already been calculated. At the 
completion of the adjustment procedure, a new point (Z,Y) k+1 
has been determined, and the following possible results must 
be considered: 

(a) . If all elements of Y k+1 are within their specified 
bounds, then f (X) is evaluated at (Z,Y) k+1 , and the procedure 
for determining the minimum continues in the normal manner. 




Figure 3-1. Adjustment of state 
variable to obtain a feasible 
point during the linear search 
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(b) . If any element of Y k+1 is not within its specified 
bounds, then (Z,Y ) k+1 is infeasible. Successive linear 
interpolation is performed between the last feasible point 
(Z,Y) k and the point (Z/Y )** 1 to determine the step length at 
which the nearest bound becomes active. Hence this step should 
conclude with a single state variable equal to one of its 
bounds and all other state variables within their specified 
bounds. Figure 3-2 shown (z 1/ y 1 ) 2 being out of bounds y x < 0) . 
Using successive false position, the bound point (z 1 ,y 1 ) 3 can 
be located. Supplementary tests are then performed to 
determine whether the local minimum lies at the bound or at 
some point before it. If the minimum lies at the bound, then 
the line search is terminated. If it lies before the bound, 
then the minimum has been bracketed, and refinement can 
be started to locate the minimum. 

(c) . If the procedure fails to converge in a reasonable amount 
of time, the step length is reduced (a = a/2), and a new trial 
point is generated. 
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Figure 3-2. Adjustment of state 
variable to locate a feasible 
point. 
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3.5. Algorithm of Generalized Reduced Gradient Method. 

According to the principle of GRG method, its algorithm 
is presented as the following: 


Step 1. Obtain the feasible initial points. 

Given a specified initial value of the search 
parameter a = a 0 , termination parameter e . 

Step 2. Choose a partition of X into state Y and decision Z 
variables such that d^/dY has nonzero determinant. 

Step 3. Calculation of the reduced gradient D : the direction 
of move for the independent variable z , by the 
following substeps: 

Step 3.1. Compute the reduced gradient D(z), given by: 

gpcf - g(Z) T - g{Y? (-^-) 


D(Zj)= 0 

if 

Zj = a 3 

and 

<h 

> 0 

D(z j )= 0 

if 

z j = bj 

and 


< 0 

D(z j )= g j 

otherwise 





if £ D(z) s termination parameter e, then a constrained 
relative minimum has been obtained. Otherwise, the algorithm 
proceeds to the next step. 
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Step 3.2. Compute D (y) , the modified reduced gradient, i.e. 
the (opposite) direction of move for the independent variable 
y. This direction may simply be D(y): 


D(Y) - - 


dW_ A 

dY 


je D(Z) 


Step 4. Compute a first value of the positive number a . 


and Compute 

z° + a 

D, 

and 

project it 

onto the 

parallelotope 

a j 5 z 3 

* 

to 

obtain z 1 . 

set : 

Zj a = a 3 

if 

Zj ° + 

a D } 

< 


z i' * 

if 

Zj ° + 

a D 3 

> b 3 



Zj 1 = Zj 0 + a D 3 otherwise 

Step 5. Compute a feasible Z 1 corresponding to a , i.e. try to 
solve, with respect to y, the system of M equations in 
M unknowns : 

V ( z 1 , Y ) =0 

This is usually done by some iterative method 
(Newton' method ) . 

Step 6. If no speedy convergence is observed, then decrease a 
for instance, = 1/2 a) and go to step 4, with 
the same D. Otherwise, let y 1 be the solution obtained 
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for iJMz 1 , y) = 0, and Z 1 the corresponding point in 
the whole n - dimensional space. 

Step 7. If f ( Z° ) < f ( Z 1 ) , then decrease a , as above, 
and go to step 4, with the same D. If or = specific 
criterion such as 10 ~ 12 , then go to step 2. Otherwise, 
at the end of step 7, we have some feasible Z 1 , which 
satisfies : 


f ( Z 1 ) < f ( Z° ) 

Step 8. We may now, either set Z° = Z 1 and begin a new 
iteration, or try to improve the last value obtained 
for a. In doing this, we return to step 4 for any new 
value tried for a, with the same D, and eventually 
terminate step 7 with some Z 1 satisfying f(Z x )< f(Z 2 ), 
and then begin a new iteration with Z° = Z 1 , go to 


step 3. 
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Figure 3-3 Flow chart for generalized reduced 
gradient (GRG) optimization method 






















CHAPTER IV 


RELIABILITY DESIGN METHOD BASED ON 
OPTIMIZATION DESIGN METHOD 


4 -1 . Introduction 

The ultimate goal in engineering design is to produce an 
optimal structure system that satisfies the performance/cost/ 
weight/volume/speed ratio objectives as well as the system or 
component reliability constraint, which is used to account for 
uncertainty existing in different failure models, loading 
conditions, material properties, and geometric parameters. To 
deal with these uncertainties, reliability technology provides 
tools for formal assessment and analysis of such 
uncertainties. However, in order to reach the optimal design 
objective, an appropriate optimizer must be used. The 
reliability design method based on optimization design method 
(also called as Probabilistic Design Optimization [PDO] ) has 
been researched by Frangopol [24 ] ; Sorensen and Thoft- 
Christensen [25]; Nicolaidis and Burdisso [26]; Maglaras and 
Nikolaidis [27]; Torng and Yang [28]; and Onwubiko et al.[29]. 

The objective function for optimal design, OBJ (X) , is 
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subject to the following constraints: 

P (g A (X) s 0U(1 - Pi), i= 1, ... , L. (4-1) 

Where X is the vector of n random variables, P(.) denotes the 
probability of the event ( . ) , g ± represents the i-th limit 
state function, ( 1- p L ) is the reliability goal for the i-th 
constraint or failure mode, and L represents the total number 
of constraints. 

In general, there are two major difficulties for the 
design problem: (1) how to solve complex problems that 

require a computation intensive program and (2) how to reduce 
the total computational effort within the design optimization 
process. To overcome these difficulties, the proposed method 
uses an advanced mean value method (AMVj [30, 31] which has 

been illustrated to be efficient for solving reliability for 
complex problems. To improve the efficiency of computation, 
the proposed method uses an approximate function to represent 
the original complex component reliability problem [32] . In 
other words, there will be only one reliability calculation in 
each design iteration. 

4-2. Optimal Structural Design Definition : 

An optimal structural system design must be insensitive 
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to uncertainties incurred from material properties, 
environmental conditions, manufacturing variations, etc . A 
smaller system variation must be achieved in order to reduce 
the possible failure [33,34]. In other words, this optimal 
structural system design must have higher reliability or lower 
probability of failure. An optimal structural system design is 
defined as a high reliability system which not only satisfies 
the performance weight / volume / cost objective but also the 
component / system reliability constraints. 

To achieve an optimal structural system design, the first 
important thing is to have a well-defined design problem. With 
consideration of design random variables, a more optimal 
structural system can be achieved; however, the optimal design 
problem setup must be redefined. In general, this new design 
optimization problem will have an objective function to be 
minimized or maximized as follows: 

Objective : F(X,Y), (4-2) 

Subject to the reliability constraints: 

P( gi (X,YUOU(l-Pi) , i=l, . . . , L. (4-3) 


Where Y is the vector of m design random variables, X is the 
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vector of n random variables, P(.) denotes the probability of 
the event ( . ) , gi represents the i - th limit state function, 
(1 - pj is the reliability goal for i-th constraint or 
failure mode, and L represents the total number of 
constraints . 

4-3. Reliability Constraint Function Definition 

With all the random variables or design random variables 
defined, different failure mechanisms - e.g., yield failure, 
fracture failure, and so on - need to be established. 
Reliability constraint functions or limit state functions are 
used to represent these failure mechanisms. These functions 
can be constructed through the response function, Z, 

Z = Z ( X, Y ) (4-4) 

where X represents the random variables and Y represents the 
design random variables. This Z function can be a simple close 
form function or a complicated model which requires the use of 
computer intensive program to model [30] . 

To calculate the reliability or probability of failure, 
a critical failure event must be defined. This failure event 
is defined when Z function value is less than or greater than 
a critical response value z D . In other words, the reliability 
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constraint function or limit state function, g, becomes: 

g ( X,Y ) = Z( X,Y ) - z Q (4-5) 

The limit state g(X)= 0 separates the variable space into 
"failure" and " safe" regions. When the equal chance constraint 
function becomes unequal, i . e . , g <0 or Z s z D the reliability 
or probability of failure, p f can be calculated as : 

P f =Prob (g (X,Y) <; 0) =Prob (Z (X,Y) - z 0 <; 0) (4-6) 

For each simple close formed g function, the reliability 
computation is straightforward. To calculate an implicitly 
defined g-function, however, the total computation becomes 
time consuming so that the selected probabilistic method must 
be efficient and reasonably accurate. 

4.4. Reliability Constraint Function Calculation: 

In general, the structural reliability analysis method is 
developed to solve a limit state function g(X). Given the 
joint probability density function, f x (X), the probability of 
failure can be formulated as: 


P f = P (g s 0) =J...J n f x ( X ) dx 


(4-7) 
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where Q is the failure region. This multiple integral is in 
general very difficult to evaluate even though there is a 
Monte Carlo solution that can provide a convenient but usually 
time-consuming solution. The first step in the current 
reliability analysis methods requires the transformation of a 
general dependent, random vector X into an independent, 
standardized normal vector u. The Rosenblatt transformation 
[16] has been suggested for this, when the joint distribution 
is available [35,36]. If only the marginal distributions and 
the covariances are known, a transformation can be made to 
generate a joint normal distribution that satisfies the given 
correlation structure. 

By transforming g(X) to g(u), the most probable point 
(MPP ) in the u-space, u*,is located. u* is the point that 
defines the minimum distance, (i, from the origin ( u = 0 
point ) to the limit state surface. This point is most 
probable because it has maximum joint probability density on 
the limit state surface, as shown in Figure 4-1. The MPP may 
be found by using optimization method or advanced mean value 
(AMV) method. Next, the g(u) or g(X) function is approximated 
by a polynomial that approximates the true function in the 
vicinity of the MPP. Once the approximate function is 
obtained, the associated failure probability can be computed. 
If the g(u) formulation is used, several analytical solutions 



Joint Probability Density 



Figure 4-1 Illustration of a most probable 
Doint [28] 
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are available for linear and quadratic functions [36] . For 
example, the first-order reliability method (FORM) estimate is 

P( gs 0 ) * <D (-3 ) (4-8) 

To compute the component structural reliability for complex 
problems that require computation intensive programs, the 
Advanced Mean Value method (AMV) is suitable because it was 
developed to search for the MPP with fewest extra g function 
calculations by comparing with the conventional mean based 
second moment method [30,31,32]. 

Let us assume that the Taylor's series expansion of 
performance function, Z, exists at the mean values. The Z 
function can be expressed as: 


Z(X) - Z(n) ♦ £ QT t - a) ♦ Hi X) 
m dX ( 

■ «o* E 
1-1 

= Z,(X) ♦ H(X) (4-9) 


where the derivatives,^, are evaluated at the mean values, p, Z : 
is a random variable representing the sum of the first-order 
terms, and H(X) represents the higher-order terms. In general, 
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the coefficients a T can be computed by numerical 
differentiation, and the minimum required number of Z-function 
calculations is (n+1) for n random variables. Since Z x is 
explicit and linear, its cdf (cumulative distribution 
function) can be computed efficiently. 

For nonlinear Z-functions, the solution based on Z x is 
only approximate. To improve accuracy, higher-order 
approximation functions can be developed. However, for 
problems involving implicit Z-functions and a large n, the 
higher-order approach might be difficult and inefficient. 

The AMV method reduces the truncation errors by replacing 
the higher-order terms H (X) by a simplified function H(Z X ) 
dependent on Z l . Ideally, the H(Z X ) function should be based 
on the exact most probable point (MPP) locus of the Z function 
to minimize the truncation error. The AMV procedure simplifies 
this approach by using the MPP of Z x . 

At each calculated MPP, probability sensitivity 
factors, a, for every defined design random variables or random 
variables are the by-product from the reliability analysis. 
These sensitivity factors, as discussed, are defined in the 
transformed standard normal space ( u - space ) : 


_ ap 

du d<t> > (F x (X)) 


( 4 - 10 ) 
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where 3 represents the safety index value, S*' 1 represents the 
inverse standard normal cdf, and F x (X) represents the cdf for 
the original random variable, X. Comparing the absolute values 
for these sensitivity factors shows their relative importance 
to the reliability solution. If all design random variables 
have uncertain means (or standard deviations), the reliability 
itself becomes a random function of these uncertain 
parameters. To measure the effect caused by these uncertain 
parameters, the probabilistic sensitivity factors, with 
respect to these uncertain design parameters and reliability 
(safety index, 3)/ can be derived as follows: 

it . it (4-11) 
du x d\ i x d\i x 


ap_ 

da x 


8p du x 

du x do x 



(4-12) 


where u x and a x represent mean and standard deviation values 
of random variable X, respectively. Since u x is a function of 
]i x and a x , 3u x / 3p x and 3u x / do x can be derived also. 

With these du x / dp x and du x / do x values evaluated, it is 
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possible to construct an approximate reliability constraint 
function, g A (y x ,a x ), as 


. °*) * Po * E 7- (c, ■ K.) 

>1 d]i xj J * 

♦ E -p- (°x t - °x w ) - ♦'‘c 1 - -p) 

jb.i 3o Xi 1 u 

■ C„. E c,n . E C,o - ♦-'(! -p) (4-13) 

>1 J k-1 


where 3o is the safety index result, p xj0 is the j-th initial 
mean value, a xk0 is the k-th initial standard deviation value, 
<1> -1 ( . ) represents the inverse normal cumulative distribution 
function (cdf ) , (1-p) is the select reliability goal, and C 0 , 
Cj , and C k are constant. 

Torng and Yang [28] have shown a safety index approximate 
function 3i (X) as : 


h*> - p. * e t®- - m. > 

>1 du J J 9 


ap 


E . 

*. 1 do 


(a - o ) 

v *i *uT 


<V 


E c ju ,• E c, 

>1 1 fc-l 


(4-14) 
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Therefore, an approximate reliability constrained function (or 
a limit state function), g A (X), can be defined as follows: 


*.(*) * 


E 

>1 


c , • E 




C k a xt - <j>'(l - p ) 


(4-15) 


or 


g A (X) * p,(X) - ♦- , (1 - p) 


(4-16) 


In order to compute efficiently a reasonably accurate system 
sensitivity, the simplest and most efficient strategy is 
accomplished first by constructing an approximate function at 
the MPP of each bottom event. By using all approximate 
functions instead of the original complex failure models, the 
probability sensitivity, with respect to mean value and 
standard deviation, can be derived as du x /d\x x and 5u x /do x t 
respectively. This sensitivity can be calculated by perturbing 
all design random variables. Total computational effort is 
reduced greatly because these perturbation analyses are 
performed based on analytical approximate functions. 

Wu, Torng, and Yang [28,31] point out AMV method is the 
best strategy for identifying the MPP for each of the failure 
models. Therefore, by using AMV, MPPs can be identified, and 
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approximate functions can be established. Once approximate 
functions are obtained, based on the objective function and 
all approximate functions, an optimization program is used to 
find the optimal design values. In this project GRG 
optimization program, which was developed by Dr. C. Onwubiko, 
is employed as an appropriate optimizer to obtain the optimal 
design values and safety index, the latter is based on the 
equation (2-2) with proper constrained functions. 

4-5 Algorithm of Reliability Design Method Based on 
Optimization Design Method 

Step 1. Define the optimal structure system or components 
requirement and construct an optimization design 
problem. 

Step 2. Use the old design point as the initial design point. 
(In general, use the mean value of design random 
variables as the initial design point. ) 

Step 3. Construct an approximation function for component 
constraint function at the initial design point 
by the following steps: 

a) . Evaluate the safety index , |3 0 , for the i-th 

reliability constraints at X c by using advanced mean 
value (AMVj method or an appropriate optimizer. 

b) . Construct an approximate constraint function for 
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component reliability constraint based on the safety 
index with respect to the main values and standard 
deviations of those design random variables. 

c) . Construct an approximate reliability constraint ( or 

limit state function ) . 

d) . Construct other approximation functions for other 

constraint functions. 

Step 4. Based on the objective function and all approximation 
functions, GRG (GRG = Generalized Reduced Gradient 
Method) .optimizer is used to find the optimal design 
values . 

Step 5. Repeat steps 2-4 until the number of iteration is 
reached or the convergence criterion is met. 




Figure. 5-2. Flow chart for reliability design 

method based on optimization method 








CHAPTER V 


DESIGN OF A GEAR TRAIN USING A RELIABILITY 
BASED OPTIMIZATION METHOD 


5-1. Introduction 

In engineering designs, the high reliability and minimum 
volume/cost/special design requirement of components / system 
is a goal pursued by engineers. A design of a gear or gear 
train is always considered an important and complex part of 
mechanical engineering design. Gears are often built into 
machines, e.g. as part of a gearbox. Smaller gears would imply 
a smaller gearbox, which leads to further savings. Tucker [37] 
says that maximizing load capacity for a given material and 
size generally results in lowest cost per horsepower 
transmitted. Willis [38] states that "weight reduction usually 
means volume reduction, which in turn lowers cost of 
materials, handling and shipping." It can be seen then that a 
good strategy is to minimize the size of the gear, not only 
because of direct saving on the gear, but also on related 
operations. Dudley [39] states, "It is often possible to 
reduce by half the length, width, and height of a gearbox by 
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simply changing from steel gears with a low hardness value to 
full-hard gear teeth. This is an 8:1 reduction in gear weight, 
which means substantial savings in material, machinery, 
storage, and shipping costs of the gearing and the housing." 
In the automobile industry smaller gears mean lighter gears; 
hence, lighter gearboxes and ultimately lighter cars. The 
search for increased efficiency (i.e., fuel economy) makes 
reducing the size of gears important regardless of initial 
cost. In the case of helicopters, reduction in size and 
weight can result in an increase in payload [40] . 

This chapter describes the minimization of the rotation 
output of the gear system, a special design requirement, using 
reliability based optimization method. 


5-2 . Model Formulation 

In designing gears, there are at least two major causes 
of failure models that are of primary concern: bending and 
contact stress. A gear train must satisfy the rotation output 
objective as well as the system reliability constraints, which 
are used to account for uncertainty existing in different 
failure models, material properties, and geometric parameters 
of gears. According to optimization design methods, the 



deterministic problem is stated as 


Minimize F v (X) = n in n(Xi/ X i+1 ) 

i = 1,2,3, .. .,M (5-1) 

subject to 


* B , J , 


°M * 0 


(5-2) 


W. 


"P i 


a i 


cMe * *> a - nb (i • 4; 




(5-3) 


N 


P i 


'O i 


Xi <; X s ; Ni * X <; N u ; 


where 

M = number of gears and pinions in system. 

Bi = gear face width (in) . ( B A = A X A ) 

P A = diametral pitch (number of tooth/in) . 

X A = pitch diameter for gears and pinions (in) . 
A = coefficient of gear face width. 

0 = the pressure angle. 
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W t = the transmitted load. 

= geometry factor of gear and pinion. 

K v = dynamic factor. 

Ppi = poisson's ratio of pinion. 
p G i = poisson's ratio of pinion. 

E Pi = modules of elasticity for pinion. 

E Gi = modules of elasticity for pinion. 
a bi = allowable bending stresses. 
o ci = allowable bending stresses. 

R pi = radius of involute on pinion. 

R Gi = radius of involute on gear. 
n in = input of rotation speed (1/min) 

N ± = number of tooth for gears and pinions (N i =D i p^ 

X,N are a column vector with n rows and the subscripts 1 and 
u represent the lower and upper bounds on X,N respectively. 

Because the probabilistic design is concerned with 
probability of failure or the reliability of system, the 
probabilistic equivalent formulation of (5-1) ; (5-2) ; (5-3) can 
be written as: 

Minimize F V (X) = n in n( X i/ x ±+i) 

M 


i — 1/2,3,..., 


(5-4) 



62 


subject to 


P [Gj, (X) £ 0]i Pi 

i = 1, 2, 3, . . . ,M (5-5) 


where 


GfX) . 


B ( J 


’b i 


(5-6) 


G«C*) - 


w. 


R 


P i 


*0 i 


0080 * B t (1 - til) (1 - nl<) 


c i 


(5-7) 


N 


P i 


"O i 


and X is a vector of n random variables and p t is the 
specified reliability level of the system. 

In terms of the principle of reliability design method 
based on optimization techniques, the formulation given in 
equations (5-4) , (5-5) , (5-6) and (5-7) were recast for 
application of reliability based optimization. They were 


Minimize F V (X) = n in n( x i/ x i+i) 
= 1, 2, 3, ... ,M 


i 


(5-8) 
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subject to 

G t (x) = Pi - rMpi) 

i = 1,2, 3, . . ,,M (5-9) 

where G L (X) are defined by equations (5-6) and (5-7) and p' 1 ( . ) 
is the inverse of the standard normal distribution function. 
Of course, it may be necessary to scale (5-9) to avoid problem 
when using nonlinear program for constrained optimization of 
the type presented in [28] [32]. 

The mean value of the pitch diameter of teeth in the 
pinions and gears is to be determined for a minimum rotation 
output of a gear train to satisfy some specified reliability 
level. It is assumed that all material properties reported are 
at their mean values. Since actual data are generally not 
available, the standard deviation, a, may be estimated by 
coefficient of variation. 


cov - 


o. 

1 



(5-10) 


where COV is the coefficient of variation, o L is the standard 
deviation, and X is the mean value of a random variable. 
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5-3 Easamgpl® Design 

To demonstrate the application of reliability based 
optimization method, we consider the following problem: 

Designing a spur gear train shown on Figure 5-1 involves 
minimizing the rotation output while satisfying the stress 
constraints. It is delivered transmitting 100 hp with a shaft 
rotating input at 2000 rmp. The material to be used is AISI 
1095. The material properties and other information are given 
in Table 5-1. 

To execute this design problem, certain assumptions are 
made. The dynamic factor k v is assumed to be 1. The J-factor 
is computed using the fitted equation given by Carrol and 
Johnson [41] . Because of the limitations of design geometric 
and undercutting, N x is assumed to be 17, and N u is assumed to 
be 100. Since the maximum contact stress occurs at the 
lowest point of single tooth contact [42], this point is 
close to the pitch point; thus, the sliding velocity is 
small . 

Therefore, the formula (5-7) is modified as follows: 


l w . ( x » 


X . x r 




1 - |1 J 1 - n 2 


a. < 0 

Ct 


( 5 - 11 ) 


N 
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Figure 5-1. Gear train for design example [43] 




Table 5-1. 


Summary of material properties and other information 
for the example problem 


Variables 

Values 

Power (pw) 

100 hp 

Rotation input 

2000 (1/min) 

Yield stress a Y i 

83 kpsi ( 572.4 MPa) 

Yield stress 

83 kpsi ( 572.4 MPa) 

Tensile stresss an 

142 kpsi ( 979.3 MPa) 

Tensile stresss a-n 

142 kpsi ( 979.3 MPa) 

Possion’s ratio \ 1 

0.30 

Pressure angle 0 

20° 

Coefficient of width X, 

0.60 

Coefficient of variation ( COV ) 

0.05 

Power efficiency r\ 

99.0 % 

Modulus of elasticity Epi ; Ep 2 

30 x 10 6 Psi ( 205 Gpa ) 

Modulus of elasticity E gJ ; 

30 x 10 6 Psi ( 205 Gpa ) 
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where in terms of Figure 5-2 


1 1 1 + 1 
R Pl R a, r Pl *** r a, 


and 




( 5 - 12 ) 


( 5 - 13 ) 


therefore 


J_ _1_ _2_ A - fo, 

R „ + *o, X 9 X ai 


( 5 - 14 ) 


After the modification of functions, then GRG 
(Generalized Reduced Gradient method) optimizer is applied to 
calculating the probability of failure, reliability and safety 
index. The probability of failure, reliability, and safety 
index are illustrated in Table 5-2, which are based on the 


calculation using GRG computer program. 
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This system can be considered to be series system. A 
series system is one in which all components are so 
interrelated that the entire system will fail even if any one 
of its components fails. Let us suppose that the components 
are independent, namely, that the performance of any one part 
does not affect the reliability of the others. Under these 
conditions, the reliability of this system is defined as 
follows : 


*,-n «> < 5 - 15 > 

M 

where R 3 is the reliability of system, R ± is the reliability 
of each component, and i is equal to 4 for this design. 

The approximate function can be constructed for the 
computation of optimization method. It is 

G ± (X) = - cjrMPi) 

i = 1,2 (5-16) 

where G ± (X) is defined by (5-6), (5-11) and p ± is defined by the 
reliability of system R s . 

After the construction of the approximate functions, an 
appropriate optimizer, GRG computer program is used in order 
to obtain the best design results for this design system. 
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Optimization functions are based on objective function (5-8) 
and constraints (5-16) . These values, rotation input n in ; 
Possion's ratio Modulus of elasticity coefficient of 
width X L ; delivered transmitting hp,are kept constant during 
the optimization process. Diametral pitch p A are assumed as 6 
(teeth/in) and 5 (teeth/in) . The results for this design 
system using reliability based optimization methods are shown 
in Tables 5-2 through 5-5 and Figures 5-3 through 5-7. All of 
calculation for this project in detail is in Appendix A. 
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Table 5-2. Results for the reliability of each component 


Situation 

Prob.of failure Pf 

Reliability R;(%) 

safety index Pi 


0.15900 x 10' 3 

99.9841 

3.599571 


0.70500 x 10' 5 

99.999295 

4.3433106 

bending case for pinion 2 
g 3 (X) function)Eq.5-2,i=2) 

0.54810 x 10‘ 2 

99.4519 

2.543948 


0.72686 x 10' 2 

99.273136 

2.443667 

system 

0.12874 x 10' 1 

98.7126265 i 

2.230000 
















Table 5-3. Results for reliability vs. the change of pitch diameter for each gear, 
rotation output ( System safety index 3 = 2.23. ) 


reliablity 

(%) 

d! , pitch 
diameter(in) 

d 2 , pitch 
diameter(in) 

d 3 , pitch 
diameter(in) 

dU , pitch 
diameter(in) 

rotation 

output(l/min) 

* 

4.03 

15.28 

6.35 

14.98 

223.71 

93.32 

- 

- 

- 

- 

- 

96.41 

3.12 

14.11 

5.24 

13.87 

166.93 

97.13 

3.28 

14.46 

5.48 

14.17 

175.73 

97.73 

3.49 

14.80 

5.75 

14.47 

187.43 

97.98 

3.60 

14.95 

5.90 

14.62 

194.33 

98.715 

4.03 

15.28 

6.35 

14.98 

223.71 

98.78 

4.08 

15.39 

6.41 

15.02 

226.15 

98.81 

4.11 

15.42 

6.45 

15.04 

228.60 

99.18 

- 


- 

- 

- 


Note: 


* indicates deterministic solution . 
- shows no feasible solution . 
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Table 5-4. Results for reliability vs. the change of the number of teeth , face width 
( System safety index p = 2.23. ) 

( Assumed the diametral pitch Pi = 6 teeth/in; P 2 =5 teeth/in) 


reliablity 

(%> 

Ni , number 
of teeth 

N 2 , number 
of teeth 

N 3 , number 
of teeth 

N 4 , number 
of teeth 

face width for 
pinion 1 ( in ) 

face width for 
pinion 2 ( in ) 

* 

24 

92 

32 

75 

2.4 

3.81 

93.32 

- 

- 

- 

- 

- 

- 

96.41 

19 

85 

26 

70 

1.87 

3.14 

97.13 

20 

87 

27 

71 

1.97 

3.29 

97.73 

21 

89 

29 

72 

2.09 

3.45 

97.98 

22 

90 

30 

73 

2.16 

3.54 

98.715 

24 

92 

32 

75 

2.40 

3.81 

98.78 

25 

92 

32 

75 

2.45 

3.85 

98.81 

25 

93 

32 

75 

2.47 

3.87 

99.18 

- 

- 

- 

- 

- 

- 


Note: * indicates deterministic solution . 

- shows no feasible solution . 


Table 5-5. Results for reliability vs. applied stresses and rotation output 
( System safety index P = 2.23. ) 


reliablity 

(%> 

bending stress 
for pinion 1 
( ksi ) 

contact stress 
for pinion 1 
(ksi) 

bending stress 
for pinion 2 
(ksi ) 

contact stress 
for pinion 2 
(ksi) 

rotation 

output 

(1/min) 

* 

10.07 

79.46 

12.17 

80.20 

223.71 

93.32 

_ 

- 

- 

- 

- 

96.41 

13.10 

78.65 

16.06 

79.00 

166.93 

97.13 

12.23 

78.90 

14.96 

79.37 

175.73 

97.73 

11.38 

79.14 

13.83 

79.72 

187.43 

97.98 

10.99 

79.24 

13.27 

79.88 

194.33 



79.46 

12.17 

80.20 

223.71 

98.78 

9.98 

79.48 

12.16 

80.22 

226.15 1 

98.81 

9.90 

79.50 

12.05 

80.23 

228.60 

99.18 

- 

- 

- 

- 

- 


Note: 


* indicates deterministic solution . 
- shows no feasible solution . 



























































































rotation output (1/min) vs. reliability (%) 



Figure 5-3. The comparison of rotation output vs. relaibility 









the changes of the number of tooth vs. reliability 



Figure 5-5, The comparison of the changes of number of tooth vs. reliability 





applied stress ( Ksi ) vs. reliability for pinion 1 



Figure 5-6. The comparison of applied stress (Ksi) vs. reliability 



applied stress ( Ksi ) vs. reliability for pinion 2 



( jsyj ) ssajjs 


Figure 5-7. The comparison of applied stress (ksi) vs. reliability (pinion 2) 
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5-4. Discussions 

To aid in the discussion of the research, the actual 
stresses in the pinion tooth are shown in Figures 5-6 and 5-7. 
The system safety index shown in Tables 5-3 through 5-5 is the 
minimum possible, based on the given information in Table 5-1 . 
By comparing the results, the rotation output of the gear 
train is increasing when the system reliability is increasing. 
This means if the higher reliability of system in design is 
selected, the heavier, larger system has to be taken. However, 
there is the limitation of a system reliability taken by a 
designer. It is impossible to obtain higher reliability of 
system more than 99.379 percent in this design. The reason is 
that reliability for an engineering system depends on the 
mean and standard deviation or Coefficient of Variation (COV) 
of design parameters. According to the principle of 
Probabilistic Design Method (PDM ) , if the higher standard 
deviation or COV in design is used, then the system will have 
the suitable reliability. Once this reliability exceeded the 
limitation of design situation, all of the design parameters 
are infeasible as shown in Figures 5-3 through 5-7 . 

From the Tables 5-3 through 5-5, the values of the design 
variables, which obtained by using Optimization Method and 
Reliability Design Method based on optimization techniques, 
are the same in the system reliability of 98.715 percent. It 
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is true that the deterministic values are the same as the 
probabilistic approach. By transforming the limit state 
functions g (X) to g(u) that had been mentioned previously, the 
most probable point (MPP) in the u-space is on the objective 
and constraint functions. Therefore, when the same safety 
index is taken, the values of design variables obtained by 
using both methods should be the same. 

Finally, we can see from the Figures 5-6 and 5-7, with 
the increase of reliability, the bending stress tends to 
decreasing while the contact stress inclines to increasing. 
This situation indicates that the failure of contact stress is 
more sensitive than the failure of bending stress. In 
conclusion, the design mainly has to be concentrated on the 
contact stress when trying to deliver high power and high 
rotation in the design of gear. 



CHAPTER VI 


CONCLUSIONS AND FUTURE RESEARCH SUGGESTIONS 


6-1. Conclusions 

The Probabilistic Design Method (PDM) is widely used in 
engineering design. PDM can be employed for stochastic design 
parameters to obtain the values of design variables under a 
specific reliability of component / system. PDM eliminates the 
deterministic design method's defect that the design variables 
must be deterministic. Also, this method makes a wider range 
of the values of design variables that can be selected by 
design engineers. However, if the design objective function is 
minimized for the volume or cost of component / system or a 
special design requirement, the values of design variables 
obtained by using PDM could not be the optimal design points 
in the design. Furthermore, this method reaches the failure of 
probability or a reliability of component / system, which is 
based on the mean and standard deviation of random design 
parameters . 

Optimization design method is a powerful tool in 
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engineering design. It emphasizes how to obtain the values of 
design variables that make the design objective function 
minimum or maximum using mathematical tools. The values of 
design variables obtained using optimization design method are 
optimal and critical; in other words, these values of design 
variables must make the system or components higher 
reliability or lower probability of failure. However, 
optimization design methods belong to a deterministic design 
method. In practice, design variables cannot be considered 
deterministic but stochastic. In addition, in the 
deterministic approach, random effects are ignored. 

There is no question that both PDM and optimization 
design method have their own disadvantages. Probabilistic 
design method was not concerned with the minimum or maximum 
design objectives but the probability of failure or the 
reliability of a system. In optimization design method, 
design variables must be deterministic. The method of 
reliability based on optimization design method eliminates 
these disadvantages with PDM and optimization design method. 
Not only that it can be applied to the situation of uncertain 
design variables in engineering design but also provide the 
wide range of reliability that designers can choose for 
engineering system or components in optimization design. 
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6-2 Suggestions for Future Research: 

Based on the results of the present study, the following 
topics may be addressed in future research: 

In this design of gear train, only both bending and 
contact stress were taken as design limit functions. Actually, 
it is complex to practically design a gear train, especially 
when the input power is more than 75 KW. Therefore, various 
design factors may be considered in future design of a gear 
train, such as thermal conditions, wear and dynamic factors, 
and so on. 

Since the safety index is defined as the minimum 
distance from the origin to the surface of the limit state 
function, minimizing the safety index using optimization 
techniques is a kind of calculated method. Therefore, there 
may be the comparison of using optimization techniques and 
NESSUS code to compute safety index and reliability of each 
component and system in future research. 
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APPENDICES 



APPENDIX A: 


A-l. The functions for calculating the safety index p 
Objective function: 


safety index { P) « ( £ x* ) 2 

/-i 


where n * number of design variables in function (slack variables are not included) 
Constrained functions: 


1) Bending function (Equation 5-2, i = 1) 

G(l)= 3.49*(1+ 0.05* x(6))- 

(l+0.05*x(5))*(1.76*x(l)*x(2)*x(3)+17.36*x(2)+6.68*x(l)) 
- (x(l) 3 x(2)*x(4))-x(7) 2 
G(2)= 1.75-x(5)-x(8) 2 
G(3)= x(6)+1.8-x(9) 2 
G(4)= x(7) 

G(5)= x(8) 

G(6)= x(9) 

where x(l)= pitch diameter of pinion 1 . x(2)= pitch diameter of gear 1 . 

x(3)= gear face width. x(4)= diametral pitch. 

x(5)= the transmitted load. x(6)= allowable bending stress. 

x(7),x(8),x(9) are slack variables. 


2) Bending function (Equation 5-2, i = 2) 


G(l)= 4.005*(1+ 0.05* x(6)>- 

x(8)*(l+0.05*x(5))*(1.76*x(l)*x(2)*x(3)+17.36*x(2)+6.68*x(l)) 
+ (x(l) 3 x(2)*x(4)*x(7))-x(9) 2 
G(2)= 1.75-x(5)-x(10) 2 
G(3)= x(6)+1.8-x(l l) 2 
G(4)=x(9) 

G(5)=x(10) 

G<6)=x(ll) 



where x(l)= pitch diameter of pinion 2. 
x(3)= gear face width. 
x(5)= the transmitted load. 
x(7)= pitch diameter of pinion 1 . 
x(9),x(10),x(l 1) are slack variables. 


x(2)= pitch diameter of gear 2. 
x(4)= diametral pitch. 
x(6)= allowable bending stress. 
x(8)= pitch diameter of gear 1 . 


3) Contact function (Equation 5-3, i = 1) 

G(l)= 0.169*(1+ 0.05* x(5)) 2 - 

(1+0.05* x(4)) * 2 . 079 * (x( 1 )+x(2)) + (x(l) 3 x(3)*x(2))-x(6) 2 
G(2)= 1.75-x(4)-x(7) 2 
G(3)= x(5)+1.32-x(8) 2 
G(4)= x(6) 

G(5)= x(7) 

G(6)= x(8) 

where x(l)= pitch diameter of pinion 1 . x(2)= pitch diameter of gear 1 . 

x(3)= gear face width. x(4)= the transmitted load. 

x(5)= allowable contact stress. x(6),x(7),x(8) are slack variables. 


4) Contact function (Equation 5-3, i = 2) 

G(l)= 0.169*(1+ 0.05* x(5)) 2 - 

(1+0.05 *x(4))*2.058*x(6)*(x(l)+x(2)) + (x(l) 3 x(3)*x(2)*x(7))-x(8) 2 
G(2)= 1.75-x(4)-x(9) 2 
G(3)= x(5)+1.32-x(10) 2 
G(4)= x(8) 

G(5)= x(9) 

G(6)=x(10) 

where x(l)= pitch diameter of pinion 2 
x(3)= gear face width. 
x(5)= allowable contact stress. 
x(7)= pitch diameter of gear 1 . 


x(2)= pitch diameter of gear 2. 
x(4)= the transmitted load. 
x(6)= pitch diameter of pinion 1. 
x(8),x(9),x(10) are slack variables. 
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A-2. Input data for calculation of the safety index p 


items 

Bending 
Function 1 

Bending 
Function 2 

Contact 
Function 1 

Contact 
Function 2 

x(l) 

-2 

-1 

-1 

-1 

x(2) 

-1 

-1 

-1 

-1 

x(3) 

1 

-1 

-1 

-1 

x(4) 

-1 

1 

1 

2 

x(5) 

1 

1 

-1.5 

-1.5 

x(6) 

-1.5 

-1.5 

1 

-1 

x(7) 

1 

-1 

1 

1 

x(8) 

1 

1 

1 

1 

x(9) 

1 

1 


1 

x(10) 


1 


1 

x(ll) 


1 



slack variables 

x(7),x(8),x(9) 

x(9),x(10),x(l 1) 

x(6),x(7),x(8) 

x(8),x(9),x(10) 





























































95 



























































A-4 The functions for calculating the optimal values of design variables 


Objective function: 


Rotation output F(X) = 2000 * 

x(2) x(4) 


Constrained functions: 

G(l)=l-10.51- (x(5)*x(l) 3 *x(2)) 

*(10.58*x(l)*x(2)+17.36*x(2)+6.68*x(l)) +2.23-P 
G(2)=l - 1 0.4*x(2)-^-(x(6)*x( 1 )*x(3) 3 *x(4)) 

*(8.815*x(3)*x(4)+17.36*x(4)+6.68*x(3) +2.23-P 
G(3)=l-588.467-x(7)*((x(l)+x(2))-(x(l) 3 *x(2))) 05 +2.23-P 
G(4)=l-585.517-x(8)*(x(2)*(x(3)+x(4)Hx(3) 3 *x(l)*x(4))) a5 +2.23-p 


where x(l)= pitch diameter of pinion 1 . x(2)= pitch diameter of gear 1 . 

x(3)= pitch diameter of pinion 2. x(4)= pitch diameter of gear 2. 

x(5)= applied bending stress for bending function 1 . 
x(6)= applied bending stress for bending function 2 . 
x(7)= applied contact stress for contact function 1 . 
x(8)= applied contact stress for contact function 2. 

P = specific safety index selected by designer. 


A-5 input data for calculation of optimal design variables 

(assumed the diametral pitch of the first pair of gears = 6 1/in; the diametral pitch 
of the second pair of gears = 5 1/in ) 


items 

x(l) 

x(2) 

x (3) 

x (4) 

x(5) 

x (6) 

x (7) 

x(8) 

initial points 

3.7 

9.7 

5 

11 

22 

25 

120 

125 

limited maximum 

16.7 

16.7 

20 

20 

83 

83 

142 

142 

limited minimum 

2.8 

2.8 

3.4 

3.4 

5 

5 

5 

5 


note: unitforx(l)tox(4): in. 
unit for x(5) to x(8): ksi. 
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A-6 output data from calculation of optimal design variables (system safety index p=2.23) 


specific safety index p 

x(l) 

x(2) 

x(3) 

x(4) 

x(5) 

x(6) 

x(7) 

x(8) 

F(x) 

P = 15 

- 

- 

- 

- 

- 

- 

* 

- 

- 

p = 1.8 

3.12 

14.1 

5.24 

13.9 

13.1 

16.1 

78.7 

79.0 

166.9 

p = 1.9 

3.28 

14.5 

5.48 

14.2 

12.2 

15.0 

78.9 

79.4 

175.7 

P =2.0 

3.49 

14.8 

5.75 

14.5 

11.4 

13.8 

79.1 

79.7 

187.4 

P =2.05 

3.60 

15.0 

5.90 

14.6 

11.0 

13.3 

79.2 

79.9 

194.3 

p = 2.23 

4.03 

15.3 

6.34 

15.0 

10.1 

12.4 

79.4 

80.1 

223.6 

P = 2.25 

4.08 

15.4 

6.41 

15.0 

10.0 

12.2 

79.5 

80.2 

226.2 

P = 2.26 

4.11 

15.4 

6.45 

15.0 

9.90 

12.1 

79.5 

80.2 

228.6 

p = 2.28 

- 

- 


- 

- 

- 

- 

- 

- 


note: - means infeasible solution 
unit for F(X): 1/min. 
unit for x(l) to x(4): in. 
unit for x(5) to x(8): ksi. 


















































































APPENDIX B: LIST OF GRG COMPUTER PROGRAM 


DECLARE SUB FTERMIN (maintest!, x!Q, n!, ne!, pvalue!) 

DECLARE SUB CONSTRAINT (x!(), cons!(), n!, ne!) 

DECLARE SUB PENALFX (x!(), conslO, fxvalue!, pvalue!, beta!, n!) 
DECLARE SUB QNEWTON (x!0, s!(), n!, tol!) 

DECLARE SUB LINESEARCH (x!(), s!0, n!, h!, tol!) 

DECLARE SUB UPDATEDFP 0 

DECLARE SUB TERMINATION 0 
DECLARE FUNCTION OBJECTIVE! (x!(), n!) 

DECLARE SUB DIRECTION 0 
DECLARE SUB UPDATEMAT 0 
DECLARE SUB PDERIVATIVE 0 

DIM SHARED x(20), xo(20), y(20), f(20), p(20), u(20), ds(20), dr(20, 20) 
DIM SHARED delx(20), delg(20), s(20), grad(20) 

COMMON SHARED ja, jb, jc, jd, jmin, xmin, ji, fbase, iter 

COMMON SHARED ql, tol, fmin, G, 1, n, m, size, cl, fcount, Igec 

COMMON SHARED mloop, fee, d, fo, h, gtest, start, hstep, update 

COMMON SHARED pvalue, ne, beta 

CLS 

PRINT 

PRINT 

INPUT "How many variables"; n 
INPUT "Total number of constraints"; ne 
DIM SHARED b(n, n), cons(ne) 

INPUT "what is the tolerance"; tol 
d = .1: fcount = 0: fee = 0: mloop = 1 
DIM SHARED xt(n), ao(n, n), gx(n) 

FOR i = 1 TO n 
PRINT "x("; i; ")": INPUT x(i) 
xo(i) = x(i): y(i) = x(i) 

NEXT i 

INPUT "enter output file please:", outfileS 
IF outfileS = "" THEN END 
OPEN outfileS FOR OUTPUT AS #2 
beta = 1 : maintest = 1 
DO UNTIL maintest = 0 

CALL QNEWTON(x(), sO, n, tol) 

CALL FTERMIN(maintest, x(), n, ne, pvalue) 
beta = 1.5 * beta 

FOR i = 1 TO n: PRINT "x("; i; ")="; x(i): NEXT i 



PRINT " pvalue="; pvalue, beta: ' INPUT o 
LOOP 

finish = TIMER 
PRINT "finish^"; finish 
maxtime = finish - sstart 

finish! = TIMER 

PRINT #2, "Program took"; maxtime; " sec "; giter; "iteration(s)" 
PRINT 

PRINT #2, "The solution after"; iter; "iteration(s) 

FOR i = 1 TO n 

PRINT #2, "x("; i; ") ="; x(i); "grad("; i; ")="; grad(i) 

NEXT i 

CALL PENALFX(xO, consO, fxvalue, pvalue, beta, n) 
fo = OBJECTIVE(xO, n) 

PRINT #2, "Objective Function="; fo, "fcount-'; fcount 
IF update = 1 THEN 

PRINT #2, "Using the DFP update and initial step of'; cstep 
ELSE 

PRINT #2, "Using the BFGS update and initial step of'; cstep 
END IF 

PRINT #2, "The solution after"; iter; "iteration(s) 

FOR i = 1 TO n 

PRINT #2, "x("; i; ") ="; x(i); "grad("; i; grad(i) 

NEXT i 

PRINT # 2 , "Objective Function="; fo, "fcount="; fcount 
IF update = 1 THEN 

PRINT #2, "Using the DFP update and initial step of'; cstep 
ELSE 

PRINT #2, "Using the BFGS update and initial step of'; cstep 
END IF 

PRINT #2, "tolerance="; tol 
END 

SUB BFGS 
DIM at(n), af(n, n) 

FOR i = 1 TOn 

delg(i) = grad(i) - gx(i) 
debc(i) = x(i) - xo(i) 

NEXT i 

'Compute first denuminator 
denuml = 0 
FOR i = 1 TOn 
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denuml = denuml + delx(i) * delg(i) 

NEXT i 

' Compute the second denominator 
denum2 = 0 
FORj = 1 TO n 
sum = 0 

FORk = 1 TOn 
sum = sum + delg(k) * ao(k, j) 

NEXT k 

denum2 = denum2 + sum * delg(j) 

NEXT j 

IF denuml o 0 AND denum2 o 0 THEN 
FOR i = 1 TO n 
FORj = 1 TO n 

b(i, j) = delx(i) * delx(j) / denuml 
NEXT j 
NEXT i 
FORj = 1 TO n 
suml = 0 
FOR k = 1 TO n 
suml = suml + ao(j, k) * delg(k) 

NEXT k 
at(j) = suml 
NEXT j 

FOR i = 1 TO n 
FORj = 1 TO n 
af(i, j) = at(i) * at(j) / denum2 
NEXT j 
NEXT i 

' Form updated matrix 
FOR i = 1 TO n 
FORj = 1 TO n 

ao(i, j) = ao(i, j) + b(i, j) - af(i, j) 

NEXT j 

NEXT i: start = 1 
ELSE 
start = 0 
END IF 
END SUB 

SUB CONSTRAINT (x(), consO, n, ne) 

REM contact function 1 

'cons(l) = .169 * (1 + .05 * x(5)) A 2 - (1 + .05 * x(4)) * 2.079 * ((x(l) + x(2)) / (x(l) A 3 * x(3) 
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*x(2)))-x(6) A 2 
’cons(2) = 1.75 - x(4) - x(7) A 2 
’cons(3) = x(5) + 1.32 - x(8) A 2 
'cons(4) = x(6) 

*cons(5) = x(7) 

'cons(6) = x(8) 

REM bending function 2 

cons(l) = 4.005 * (1 + .05 * x(6)) - x(8) * (1 + .05 * x(5)) * (1.76 * x(l) * x(2) * x(3) + 17.36 * 

x(2) + 6.68 * x(l» / (x(l) A 3 * x(2) * x(4) * x(7)) - x(9) A 2 

cons(2) = 1.75 - x(5) - x(10) A 2 

cons(3) = x(6) + 1.8 - x(l 1) A 2 

cons(4) = x(9) 

cons(5) = x(10) 

cons(6) = x(l 1) 

REM bending function 1 

'cons(l) = 3.49 * (1 + .05 * x(6)) - (1 + .05 * x(5)) * (1.76 * x(l) * x(2) * x(3) + 17.36 * x(2) + 

6.68 * x(l)) / (x(l) A 3 * x(2) * x(4)) - x(7) A 2 

'cons(2) = 1.75 - x(5) - x(8) A 2 

*cons(3) = x(6) + 1.8 - x(9) A 2 

’cons(4) = x(7) 

'cons(5) = x(8) 

’cons(6) = x(9) 

REM contact function 2 

•cons(l) = .169 * (1 + .05 * x(5)) A 2 - (1 + .05 * x(4)) * 2.06 * x(6) * ((x(l) + x(2)) / (x(l) A 3 * 

x(3) * x(2) * x(7))) - x(8) A 2 

'cons(2> = 1 .75 - x(4) - x(9) A 2 

’cons(3) = x(5) + 1.32 - x(10) A 2 

’cons(4) = x(8) 

’cons(5) = x(9) 

'cons(6) = x(10) 

END SUB 

SUB DIRECTION 
IF start = 0 THEN 
'Set the identity matrix 
FOR i = 1 TOn 
FOR j = 1 TO n 
IF i = j THEN 
ao(i,j)= 1 
ELSE 
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ao(i, j) = 0 
END IF 
NEXT j 
NEXT i 
END IF 

'Identify point at which A is set to identiy matrix 
lgec = gee 

' Save the initial point and the gradient at this point 
FORj = 1 TO n 
xo(j) = xQ 
gx(j) = gradO) 

NEXT j 

' Form the step and product of step and gradient 
cl = 0: snorm = 0 
FOR i = 1 TO n 
sum = 0 

FORj = 1 TO n 

sum = sum - ao(i, j) * grad(j) 

NEXT j 

s(i) = sum: cl = cl - sum * grad(i) 
snorm = snorm + s(i) A 2 
NEXT i 

'Check if normalization is necessary 
IF snorm >100 THEN 
snorm = snorm A .5 
FOR i = 1 TO n 
s(i) = s(i) / snorm 
NEXT i 
END IF 
END SUB 

SUB FTERMIN (maintest, x(), n, ne, pvalue) 

IF pvalue <= tol OR beta > 1E+20 THEN 
maintest = 0 
ELSE 

CALL PDERJVATIVE 
gradvalue = 0 
FOR i = 1 TOn 

gradvalue = gradvalue + grad(i) A 2 
NEXT i 

gradvalue = gradvalue A .5 
IF gradvalue <= tol THEN 
maintest = 0 
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END IF 
END IF 
END SUB 

SUB LINESEARCH (xO, s(), n, h, tol) 

DIMxl(n) 
toll = .0001 

CALL PENALFX(xO, consO, fxvalue, pvalue, beta, n) 
fl = fxvalue 
FORt = 1 TOn 
xl(t) = x(t) 

NEXT t 
FORt = 1 TOn 
x(t) = xl(t) + h * s(t) 

NEXT t 

CALL PENALFX(xO, consO, fxvalue, pvalue, beta, n) 
ff= fxvalue 
IFff<fl THEN 
routel = 1 

DO UNTIL routel = 0 
f2 = ff: 

FORt = 1 TOn 

x(t) = xl(t) + 2 * h * s(t): 

NEXT t 

CALL PENALFX(x(), consO, devalue, pvalue, beta, n) 
ff = fxvalue 
IF ff > f2 THEN 
routel = 0: f3 = ff 
ELSE 
h = 2 * h 
END IF 
LOOP 
ELSE 

routel = 1 

DO UNTIL routel = 0 
f3 = ff: 

FORt = 1 TOn 
x(t) = xl(t) + .5 * h * s(t): 

NEXT t 

CALL PENALFX(x0, consO, fxvalue, pvalue, beta, n) 
ff = fxvalue 
IFff<fl THEN 
routel =0:f2 = ff:h = h/2 



ELSE 

h = -.5 * h: ' Note, you may change sign to "+" 

END IF 

IF ABS(h) <= IE- 14 THEN 
route 1 = 0 
END IF 
LOOP 
END IF 

IF h > IE-14 THEN 

d = .5 * h * (4 * £2 - 3 * fl - G) / (2 * £2 - fl - f3) 
a = 0:b = h:c = 2*h 
testl = 1 

DO UNTIL testl = 0 
FORt= 1 TOn 
x(t) = xl(t) + d * s(t): 

NEXT t 

CALL PENALFX(x(), consO, fxvalue, pvalue, beta, n) 
f4 = fxvalue 
'Check convergence 

IF ABS(f2 - f4) <= toll OR ABS(b - d) <= toll THEN 
IF f4 < £2 THEN 
alpha = d 
ELSE 
alpha = b 
END IF 

FORt= 1 TOn 
x(t) = xl(t) + alpha * s(t) 

PRINT "x("; t; x(t) 

NEXT t 

CALL PENALFX(x(), consO, fxvalue, pvalue, beta, n) 
fopt = fxvalue 
PRINT "fopt="; fopt 
testl = 0 
ELSE 

' Check that bracket is not lost in discarding the max. pt 
IF a <= d AND d <= b THEN 
IF fl >= f4 AND f4 <= £2 THEN 
c = b: 0 = £2: b = d: f2 = f4 
ELSE 

a = d: fl = f4 
END IF 
ELSE 

IF £2 >= f4 AND f4 <= f3 THEN 
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a = b: fl = f2: b = d: f2 = f4 
ELSE 

c = d: f3 = f4 
END IF 
END IF 
END IF 

num = (b A 2-c A 2)*fl+(c A 2-a A 2)*£2 + (a A 2-b A 2)*f3 
den = (b - c) * fl + (c - a) * £2 + (a - b) * D 
IF den = 0 THEN 
d = 0 
ELSE 

d = .5 * num / den 
END IF 
LOOP 
ELSE 

FOR t = 1 TOn 
x(t) = xl(t) 

NEXT t 
END IF 
END SUB 

FUNCTION OBJECTIVE (x(X n) 

OBJECTIVE = (x(l) A 2 + x(2) A 2 + x(3) A 2 + x(4) A 2 + x(5) A 2 + x(6) A 2 + x(7) A 2 + x(8) A 
2) A .5 

END FUNCTION 

SUB PDERIVATIVE 
DIMzl(lOO) 

1 Subroutine computes the gradient 
delta = .001: status = 1 
FOR j = 1 TO n 
zlG) = x(j) 
x(j) = zl(j) + delta 

CALL PENALFX(xO, consO, fxvalue, pvalue, beta, n) 
yl = fxvalue 
x G) = zlG)- delta 

CALL PENALFX(xO, consO, fxvalue, pvalue, beta, n) 
y2 = ficvalue 

grad(j) = (yl - y2) / (2 * delta) 

x G) = zlG) 

NEXT j 
END SUB 


SUB PENALFX (x(), consO, fxvalue, pvalue, beta, n) 

CALL CONSTRAINT(xO, consO, ne, n) 
sumc = 0 
FOR i = 1 TO ne 
sumc = sumc + cons(i) A 2 
NEXT i 

pvalue = sumc A .5: 

' Handle the bounds on variables 
fxvalue = OBJECTIVE(xO, n) 
fxvalue = fxvalue + beta * (sumc) 

END SUB 

SUB QNEWTON (x(), s(), n, tol) 
status = 1 : pass = 0: iter = 0 
hstep = 1 : mloop = 1 : cstep = hstep 
Evaluate gradient at current point x() 
CALLPDERIVATIVE: gee = gee +1 
'Start iteration process 
start = 0 
sstart = TIMER 

CALL PENALFX(x(), consO, fxvalue, pvalue, beta, n) 
fo = fxvalue 
fee = fee + 1 
DO UNTIL mloop = 0 
CALL DIRECTION 
1 Test for search direction 

CALL PENALFX(xO, consO, fxvalue, pvalue, beta, n) 
fbase = fxvalue 
' Perform a line search 


h = hstep 

CALL LINESEARCH(xO, s0, n, h, tol) 

CALL PENALFX(xO, consO, fxvalue, pvalue, beta, n) 
fmin = fxvalue 
iter = iter + 1 

CALL PDERIVATIVE 

CALL TERMINATION: 'Check for convergence 
IF mloop o 0 THEN 
IF fmin >= fbase THEN 

start = 0: ' Restart the search using old direction 
but reduce step size 
hstep = h / 2 

IF h <= .000001 * cstep THEN 



mloop = 0 
END IF 
ELSE 

CALL UPDATEDFP 
’Set a new search direction 
hstep = cstep 
END IF 
END IF 
LOOP 
END SUB 

SUB TERMINATION 
xtest = 0: ftest = 0: gtest = 0: gtol = .000001 
ftol = .000001: xtol = .000001 
FOR i = 1 TOn 

'IF ABS(xo(i) - x(i)) <= xtol THEN 
xtest = xtest + ABS(xo(i) - x(i)) 

’ END IF 

gradvalue = gradvalue + (grad(i)) A 2 
NEXT i: gtest = (gradvalue) A .5 
IF fbase o 0 THEN 
ftest = (fmin - fbase) 

END IF 

IF gtest <= gtol THEN 
mloop = 0 
END IF 
END SUB 

SUB UPDATEDFP 
DIM at(n), af(n, n) 

FOR i = 1 TOn 

delg(i) = grad(i) - gx(i) 
delx(i) = x(i) - xo(i) 

NEXT i 

'Compute first denuminator 
denuml = 0 
FOR i = 1 TO n 

denuml = denuml + delx(i) * delg(i) 

NEXT i 

' Compute the second denominator 
denum2 = 0 
FOR j = 1 TO n 
sum = 0 



FORk= 1 TOn 
sum = sum + delg(k) * ao(k, j) 

NEXT k 

denum2 = denum2 + sum * delg(j) 
NEXT j 

IF denuml o 0 AND denum2 o 0 THEN 
FOR i = 1 TOn 
FORj = 1 TO n 

b(i, j) = delx(i) * delx(j) / denuml 
NEXT j 
NEXT i 
FORj = 1 TO n 
suml = 0 
FOR k = 1 TO n 
suml = suml + ao(j, k) * delg(k) 
NEXT k 
at(j) = suml 
NEXT j 

FOR i = 1 TOn 
FORj = 1 TOn 
af(i, j) = at(i) * at(j) / denum2 
NEXT j 
NEXT i 

1 Form updated matrix 

FOR i = 1 TO n 
FORj = 1 TO n 

ao(i, j) = ao(i, j) + b(i, j) - af(i, j) 
NEXT j 

NEXT i: start = 1 
ELSE 
start = 0 
END IF 
END SUB 

SUB UPDATEMAT 
DIM AA(n, n), af(n, n) 

FOR i = 1 TO n 

delg(i) = grad(i) - gx(i) 
delx(i) = x(i) - xo(i) 

NEXT i 
denum = 0 
FOR i = 1 TO n 
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denum = denum + delx(i) * delg(i) 

NEXT i 

IF denum o 0 THEN 
FOR i = 1 TO n 
FOR j = 1 TO n 
IF i=j THEN 

b(i, j) = 1 - (delx(i) * delg(j)) / denum 
ELSE 

b(i, j) = -(delx(i) * delgQ) / denum 
END IF 
NEXT j 
NEXT i 
FOR i = 1 TO n 
FOR j = I TO n 
sum = 0 

FOR k = 1 TO n 
sum = sum + b(i, k) * ao(k, j) 

NEXT k 
af(i, j) = sum 
NEXT j 
NEXT i 

FOR i = 1 TOn 
FORj = 1 TO n 
sum = 0 

FOR k = 1 TO n 

sum = sum + af(i, k) * b(k, j) 

NEXT k 

AA(i, j) = sum + (delx(i) * delxQ) / denum 
ao(i, j) = AA(i, j) 

NEXT j 
NEXT i 
FOR i = 1 TOn 
FORj = 1 TO n 
NEXT j 

NEXT i: start = 1 
ELSE 

start = 0 
END IF 
END SUB 



